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1  INTRODUCTION 


1.1  Motivation 

The  aerodynamic  characteristics  (e.g.,  lift,  drag  and  moment)  of  an  airfoil  in  unsteady 
motion  are  significantly  affected  by  the  viscous  boundary  layer.  Under  certain  condi¬ 
tions,  the  boundary  layer  on  an  airfoil  in  pitching  motion  separates  from  the  airfoil 
surface,  forming  a  large  concentrated  region  of  vorticity  known  as  the  dynamic  stall 
vortex.  This  phenomenon  is  characterized  by  dramatic  changes  in  the  aerodynamic  per¬ 
formance  (e.g.,  a  rapid  change  in  the  moment)  of  the  airfoil,  and  is  of  significant  interest 
to  rotorcraft,  for  example.  However,  the  complex  unsteady  effects  have  not  yet  been 
completely  understood,  and  there  is  a  need  for  more  fundamental  research  in  this  field. 

Boundary  layer  separation  is  an  integral  part  of  the  dynamic  stall.  Qualitatively, 
boundary  layer  separation  is  the  breakdown  of  the  boundary  layer  model  which  divides 
the  flow  into  two  weakly  interacting  regions  (i.e.,  an  irrotational  flow  occupying  most 
of  the  fluid  volume  and  a  thin  viscous  region  adjacent  to  the  solid  boundary).  For  2-D 
steady  flows,  the  criterion  for  boundary  layer  separation  is  well  known,  and  indicated  by 
the  appearance  of  zero  shear  stress  at  the  surface.  For  a  pitching  airfoil,  the  breakdown  of 
the  boundary  layer  (separation)  is  preceded  by  the  formation  of  a  large  scale  recirculating 
region^  above  the  airfoil  surface.  The  recirculating  region  grows  in  size  in  both  transverse 
and  normal  directions  and  at  a  later  stage  detaches  from  the  airfoil  surface  (separation) 
to  give  rise  to  the  phenomenon  of  dynamic  stall. 

The  focus  of  the  present  research  is  the  understanding  of  the  initial  (“incipient”) 
stages  of  boundary  layer  separation  over  a  pitching  airfoil,  with  particular  emphasis  on 
the  details  of  the  separation  phenomenon  near  the  leading  edge.  Improved  understanding 
of  the  incipient  stages  of  boundary  layer  separation  may  lead  to  methods  for  modification 

recirculating  region  can  be  defined  as  a  fiow  possessing  vorticity  with  closed  streamlines. 
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or  control  of  the  separation  process. 


1.2  Literature  Survey 

A  number  of  recent  research  activities  have  focused  on  the  study  of  the  processes  lead¬ 
ing  to  the  dynamic  stall  and  the  dynamic  stall  phenomenon  itself.  Experimental  studies 
have  continued  to  provide  very  important  insights  into  the  physics  of  the  dynamic  stall 
process.  McCroskey  et  al  [24]  studied  incompressible  boundary  layer  separation  using  oil 
smoke  visualization  and  described  the  three  different  types  of  boundary  layer  separation 
observed  in  the  case  of  an  oscillating  airfoil.  Acharya  and  Metwally  [1]  categorized  and 
quantified  the  sources  of  vorticity  and  identified  the  key  mechanisms  in  the  initiation, 
development,  growth  and  movement  of  the  dynamic  stall  vortex.  Chandrasekhara  and 
Ahmed  [8]  obtained  the  instantaneous  velocity  measurements  over  an  oscillating  airfoil 
in  a  compressible  medium  which  showed  the  formation  of  the  separation  bubble  over 
the  airfoil  that  persists  until  angles  close  to  when  the  dynamic  stall  vortex  forms  and 
convects.  Carr  et  al  [5]  studied  dynamic  stall  using  real  time  interferometry  and  the 
measurements  of  the  flow  near  the  leading  edge  of  an  oscillating  airfoil  offered  the  de¬ 
tailed  experimental  quantification  of  the  locally  compressible  flow  field  that  surrounds 
an  oscillating  airfoil  at  moderate  subsonic  Mach  numbers  and  also  revealed  significant 
characteristics  of  the  complex,  and  rapidly  varying  locally  supersonic  flow.  Karim  [18] 
experimentally  studied  the  evolution  of  the  dynamic  stall  vortex  in  the  vicinity  of  lead¬ 
ing  edge  of  a  2-D  pitching  airfoil  using  smoke-wires  and  studied  the  pitch  rate  effects  on 
the  growth  of  the  dynamic  stall  vortex.  He  also  investigated  the  control  of  dynamic  stall 
using  suction.  Chandrasekhara  et  al  [9]  used  real-time  point  diffraction  interferometry 
to  study  the  compressible  dynamic  stall  over  a  airfoil  pitching  at  a  constant  rate  and 
observed  the  appearance  of  small  multiple  shocks  near  the  leading  edge  above  the  shear 
layer  for  free  stream  Mach  numbers  above  0.4.  Crisler  et  al  [11]  investigated  the  un¬ 
steady  flow  over  a  airfoil  pitching  at  a  constant  rate  using  a  Particle  Image  Velocimetry 
(PIV)  system  and  identified  the  role  of  absolute  instability  of  a  separating  shear  layer 
in  the  dynamic  stall  process. 

Computational  studies  have  long  been  used  to  improve  the  understanding  of  the 
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unsteady  flow  behavior.  Mehta  [26]  computed  the  laminar  flow  past  an  oscillating  airfoil 
at  Reynolds  numbers  (based  on  the  airfoil  chord)  of  5  x  10^  and  10^  to  gain  insight  into 
the  mechanism  of  dynamic  stall.  He  qualitatively  compared  his  computed  instantaneous 
streamlines  with  the  trajectories  of  air  bubbles^  observed  in  water  tunnel  experiments. 
Shih  et  al  [37]  observed  that  the  leading  edge  boundary  layer  separation  leads  to  the 
formation  of  a  vortical  structure  which  dominates  the  aerodynamic  performance.  Ghia 
et  al  [15]  and  Yang  et  al  [58]  analyzed  the  effect  of  modulated  suction  or  injection  in 
delaying  the  onset  of  dynamic  stall  over  pitching  NAG  A  airfoils  and  also  successfully 
compared  their  computational  results  with  the  previous  computations  of  Mehta  [26] 
and  the  experimental  results  for  Rcc  =  10^  and  Rcc  =  4.5  x  10^  .  They  monitored  the 
separation  and  through  numerical  experimentation  devised  suitable  leading-edge  suction 
and  equal  mass  of  trailing  edge  injection  (necessary  to  maintain  zero  net  normal  mass 
addition  and  to  keep  the  boundary  condition  at  inflnity  unchanged,  due  to  the  use  of 
streamfunction-vorticity  formulation  and  constant  vorticity  boundary  condition  around 
the  airfoil)  to  manage  the  evolution  of  the  dynamic  stall  vortex,  and  depicted  the  role 
of  separation  in  the  formation  of  the  stall  vortex.  Visbal  [52]  presented  investigations 
of  flow  control  by  boundary  layer  suction  and  a  moving  wall.  He  also  investigated 
the  effect  of  compressibility  on  dynamic  stall  and  determined  that  an  increase  in  Mach 
number  reduces  the  stall  delay  [51].  Currier  and  Fung  [12]  found  that  with  increasing 
unsteadiness  the  delay  between  static  stall  angle  and  dynamics  stall  onset  decreases. 
Shida  et  al  [36]  computed  the  flowfield  around  an  oscillating  NACA-0012  airfoil  by 
solving  the  two-dimensional  compressible  Navier-Stokes  equations  employing  a  block 
pentadiagonal  matrix  scheme  and  were  able  to  capture  the  lift  stall  and  lift  restoration 
in  the  downstroke  of  the  airfoil  motion. 

Analytical  studies  have  been  very  helpful  in  the  fundamental  understanding  of  un¬ 
steady  separation  and  linking  the  physics  with  the  experimental  and  numerical  obser¬ 
vations.  Doligalski  et  al  [13]  have  reviewed  some  of  the  important  studies  conducted  on 
vortex  interactions  and  separation  including  dynamic  stall.  Smith  [38,  39,  41]  described 
the  instability  of  a  leading  edge  separation  bubble  and  finite  time  breakup  of  the  bound¬ 
ary  layer.  He  found  that  initially  unsteady  developments  take  place  over  a  relatively 

^Air  bubble  trajectories  are  not  instantaneous  streamlines  in  an  unsteady  flow. 
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slow  time  scale  but  then  the  corresponding  solution  breaks  down  with  a  singularity, 
forcing  a  switch  to  a  faster  and  more  nonlinear  process.  Peridier  et  al  [28,  29]  examined 
the  interaction  of  a  vortex  with  a  boundary  layer.  Some  additional  research  has  also 
focused  on  some  geometrically  simpler  configurations  such  as  a  circular  cylinder  set  into 
motion  impulsively  [47,  40].  Improved  understanding  of  the  temporal  and  spatial  scales 
associated  with  the  dynamic  stall  process  have  been  obtained  from  these  studies. 

Comprehensive  reviews  of  advances  in  the  field  of  computational  and  experimental 
studies  of  dynamic  stall  have  been  presented  by  Carr  [4]  and  Carr  and  McCroskey  [6]. 

1.3  Present  Research 

The  main  objective  of  the  present  paper  is  improved  understanding  of  the  effects  of 
compressibility,  pitch  rate  and  Reynolds  number  on  the  initial  stages  of  unsteady  leading 
edge  boundary  layer  separation  for  a  pitching  airfoil.  An  NACA-0012  airfoil  has  been 
selected,  consistent  with  previous  computational  and  experimental  studies.  Figure  1 
shows  schematically  the  typical  stages  in  the  pitch  up  motion  of  an  airfoil,  where  the 
instantaneous  streamlines  are  displayed  in  a  reference  frame  attached  to  the  airfoil.  At 
low  angles  of  attack,  the  flow  is  attached  to  the  airfoil  surface  and  there  is  no  reversed 
flow  present  in  the  flowfield.  With  increase  in  the  angle  of  attack,  a  reverse  flow  region 
forms  over  the  airfoil  top  surface  and  finally  extends  up  to  the  leading  edge.  A  leading 
edge  vortex  forms  over  the  airfoil  which  eventually  moves  away  from  the  airfoil  surface 
to  give  rise  to  dynamic  stall.  The  present  research  deals  with  the  understanding  of  the 
stages  from  (b)  (attached  flow)  through  (e)  (dynamic  stall),  but  not  including  dynamic 
stall  itself. 

The  three-dimensional  parameter  space  of  Mach  number,  pitch  rate  and  Reynolds 
number  has  been  investigated  (where  the  Reynolds  number  is  based  on  the  chord  length 
c  and  the  non-dimensional  pitch  rate  is  0+  =  Uc/Uoo,  where  Q,  is  the  pitch  rate  in 
rad/s  and  Uoo  is  the  freestream  velocity).  The  seven  points  studied  in  the  three- 
dimensional  space  are  indicated  by  solid  circles  in  figure  2.  Some  of  the  previous  lam¬ 
inar  Navier-Stokes  simulations  of  a  pitching  airfoil  due  to  other  researchers  have  also 
been  indicated  on  the  figure.  Computations  have  been  performed  using  two  different 
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(C) 

reversed  flow 


attached  flow 


leading  edge  vortex  forms 


Figure  1:  Stages  in  the  pitching  of  an  airfoil 

numerical  algorithms.  The  first  algorithm,  denoted  the  structured  grid  algorithm,  is 
an  approximate-factorization  implementation  of  Beam-Warming’s  method  [2]  using  a 
structured,  boundary-fitted  grid  system.  The  second  algorithm  [20],  denoted  the  un¬ 
structured  grid  algorithm,  employs  an  unstructured  grid  of  triangles  and  utilizes  the 
flux- differencing  splitting  method  of  Roe  [34]  for  the  inviscid  fluxes  and  a  discrete  rep¬ 
resentation  of  Gauss’  Theorem  for  the  viscous  fluxes  and  heat  transfer.  The  structured 
grid  algorithm  is  implicit,  while  the  unstructured  grid  algorithm  (as  employed  in  this 
study)  is  explicit.  Both  algorithms  are  second  order  accurate  in  space  and  time.  The 
structured  grid  solver  has  been  employed  to  compute  all  the  seven  cases.  The  unstruc- 
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tured  grid  solver  has  been  employed  to  compute  one  of  the  cases  in  order  to  provide 
an  independent  examination  of  the  accuracy  of  the  computations.  The  airfoil  is  pitched 
about  the  quarter  chord  axis.  The  flow  conditions  for  the  seven  computed  cases  were 
chosen  on  the  basis  of  simplicity  and  feasibility.  The  Reynolds  numbers  were  selected 
to  ensure  laminar  flow.  The  Mach  numbers  were  chosen  to  obtain  subsonic  flow  in  one 
case  and  regions  of  supersonic  flow  in  the  other  case.  The  studies  were  conducted  with 
special  emphasis  on  understanding  the  leading  edge  separation,  of  course.  The  boundary 
layer  separates  near  the  trailing  edge  very  early  in  the  pitch  up  motion,  but  for  most 
of  the  cases  it  is  the  separation  of  the  boundary  layer  near  the  leading  edge  which  is 
responsible  for  the  eventual  dynamic  stall  process. 

(0.6) 

□  Visbal 
O  Ghia  et  al. 

O  Shih  et  al. 

•  present  study 


Figure  2:  Points  investigated  in  the  three-dimensional  parameter  space  (pitch  rate  value 
is  indicated  by  the  number  in  the  bracket) 

The  early  study  of  Ghosh  Choudhuri  et  al  [16]  at  Rcc  =  10“^,  Moo  =  0.2,  and  non- 
dimensional  pitch  rate  =  0.2  reported  the  appearance  of  the  primary  recirculating 
region  (a  flow  structure  possessing  vorticity  with  closed  streamlines)  to  be  a  very  impor¬ 
tant  development  leading  to  the  separation  process.  The  primary  recirculating  region 
eventually  separates  from  the  airfoil  surface  to  give  rise  to  the  dynamic  stall  vortex. 
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Therefore,  the  appearance  of  the  primary  recirculating  region  has  been  closely  moni¬ 
tored.  A  linear  stability  analysis  has  also  been  conducted  to  determine  whether  the 
appearance  of  the  primary  recirculating  region  is  related  to  the  instability  of  the  flow- 
field. 

This  study  differs  from  previous  computational  studies  in  the  following  aspects  : 

•  This  is  a  comprehensive  study  of  an  expanded  three-dimensional  parametric  space 
of  Mach  number,  pitch  rate  and  Reynolds  number. 

•  The  initial  stages  of  the  development  of  the  boundary  layer  leading  to  the  devel¬ 
opment  of  the  dynamic  stall  vortex  has  been  studied. 

•  A  linear  stability  analysis  has  been  performed  to  understand  the  reason  for  the 
formation  of  the  primary  recirculating  region. 

The  governing  equations  for  the  numerical  simulations  utilizing  structured  and  un¬ 
structured  grid  algorithms  have  been  presented  in  sections  2  and  4,  respectively.  The 
algorithms  have  been  described  in  detail  in  sections  3  and  5.  These  sections  also  contain 
details  about  the  boundary  conditions  used  for  the  numerical  simulations.  Sections  6 
and  7  discuss  the  code  validation  study  for  a  number  of  cases  through  comparison  with 
analytical  or  previous  computations.  Section  8  specifies  the  problem  itself.  Section  9 
consists  of  the  results  for  the  present  study.  Section  10  provides  a  summary  of  the 
important  results  and  future  work. 
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2  GOVERNING  EQUATION 
(STRUCTURED  GRID) 


The  governing  equations  employed  for  the  numerical  simulation  of  unsteady  flow 
past  an  airfoil  utilizing  a  structured  grid  are  presented  in  this  chapter.  The  equations  are 
the  two-dimensional,  unsteady,  compressible,  laminar  Navier-Stokes  equations  written 
in  strong  conservation  form  for  general  curvilinear  coordinates. 

2.1  Equations  in  Cartesian  Coordinates 

The  unsteady,  compressible,  two-dimensional  Navier-Stokes  equations  in  strong  con¬ 
servation  form  can  be  expressed  in  the  Cartesian  coordinates  in  the  following  form  : 


m  dF  ^  M  ^ 
dt  ^  dx  dy  dx  dy 
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The  density,  p,  static  pressure,  p,  and  the  absolute  temperature,  T,  obey  the  equation 
of  state  : 

1 


p  =  pRT  =  p(7  -  l){e  -  +  u^)} 


(3) 
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where  R  is  the  Universal  gas  constant,  e  is  the  total  energy  per  unit  mass,  A  is  the  second 
coefficient  of  viscosity  (A  =  — 2/3/x),  and  7  is  the  ratio  of  specific  heats. 

The  dynamic  molecular  viscosity  is  assumed  to  satisfy  Sutherland’s  relation  [56]: 


if 

/i 


T  To  +  Si 


(4) 


T  +  Si 

where  is  the  Sutherland’s  reference  temperature  (198°i?  for  air),  and  fio  =  y-{To)- 
The  molecular  Prandtl  number,  Pr  (0.73  for  air),  and  the  specific  heat  at  constant 
pressure,  Cp,  are  assumed  constant. 


2.2  Equations  in  General  Curvilinear  Coordinates 


The  aerodynamic  flow  over  an  airfoil  is  characterized  by  a  relatively  complex  con¬ 
figuration.  It  is  very  difficult  to  describe  such  configurations  in  the  standard  cartesian 
or  polar  coordinates.  Writing  the  governing  equations  in  terms  of  general  “boundary- 
conforming”  curvilinear  coordinates  helps  in  simplifying  the  numerical  simulation  of 
aerodynamic  flows.  The  physical  boundaries  of  the  flow  are  mapped  into  constant  trans¬ 
formed  coordinate  lines,  and  this  eliminates  the  inaccuracies  and  the  need  for  interpola¬ 
tion  in  the  numerical  implementation  of  the  boundary  conditions.  It  also  allows  a  more 
efficient  resolution  of  the  features  of  the  flow  when  high  gradients  are  associated. 

A  time  dependent  non-singular  transformation  is  introduced  :  {x,y)  (^,»/),  where 

(  =  ^(x,y,t),  and  r]  =  T]{x,y,t).  With  this  transformation  the  Navier- Stokes  equations 
(Eqn.  (1))  can  be  written  as  [43]: 
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The  Jacobian  (J)  of  the  transformation  from  {x,y)  to  (^,?7),  and  the  contravariant 
velocities,  U  and  V,  are  defined  as 
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The  coefficients  6,,  Cj,  di  (i=l,..,4;  j=l,...,5)  appearing  in  the  viscous  term  of  eqn.  (5) 
are  defined  as  follows 
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The  transformation  metrics  ^t,  7/t,  7/3.,  and  rjy  appearing  in  the  above  equations 

can  be  defined  in  terms  of  x^,  Xr,,  y^,  and  y^,  which  can  be  computed  numerically  by 
applying  a  finite-difference  formula  to  the  body-conforming  grid. 
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The  Beam- Warming  algorithm  has  been  applied  to  the  set  of  equations  presented  in 
this  section  and  the  details  of  the  algorithm  axe  presented  in  the  next  section. 
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3  NUMERICAL  ALGORITHM 
(STRUCTURED  GRID) 


Unsteady  aerodynamic  flow  problems  require  the  temporal  accuracy  along  with 
the  spatial  accuracy  in  the  solution  of  the  unsteady  Navier-Stokes  equations.  Numerical 
integration  in  time  is  performed  to  solve  the  governing  equations  and  the  transients  are 
of  prime  importance.  It  is  also  desirable  to  employ  a  time-marching,  finite  difference 
algorithm  that  is  unconditionally  stable  and  therefore  allows  the  use  of  a  large  time 
increment,  to  solve  the  steady  part  of  the  problem. 

Beam  and  Warming  [2]  have  developed  a  fully  impbcit,  unconditionally  stable,  non¬ 
iterative  finite-difference  scheme  for  the  solution  of  hyperbolic  and  mixed  hyperbolic- 
parabolic  systems  of  partial  difference  equations.  The  Beam- Warming  scheme  em¬ 
ploys  a  time-linearization  procedure  which  eliminates  the  need  for  an  iterative  method 
at  each  time  step  and  utilizes  the  ADI  technique,  which  results  in  the  solution  of 
block-tridiagonal  linear  systems  for  each  coordinate  direction.  It  utilizes  a  structured 
boundary-fitted  grid.  The  Beam-warming  algorithm  in  its  “delta”  formulation  with 
Trapezoidal  implicit  time  differencing  and  centered  spatial  approximations  is  second  or¬ 
der  accurate  in  time  and  space.  This  section  presents  the  numerical  implementation 
of  Beam- Warming  algorithm.  The  boundary  conditions  and  grid  generation  are  also 
discussed. 


3.1  Beam- Warming  Algorithm 


The  governing  equations  are  solved  in  an  inertial  frame  of  reference.  Trapezoidal 
temporal  discretization  of  the  Navier-Stokes  equations  (eqn.  (5))  gives 
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where  n  indicates  the  temporal  index  and  At  is  the  time  increment.  Therefore,  = 
q{nAt)  =  q{t).  Also, 

A5”  =  9”+^  -  g";  AjP”  =  -  F”;  AG"  =  G"+^  -  G";  AV^"  =  -  F" 

AFj”  =  AH^i"  =  AW^s"  =  (18) 


The  delta-form  of  Beam- Warming  algorithm  results  in  several  advantages,  including 
a  steady  state  solution  independent  of  At  and  a  more  direct  derivation  of  the  factored 
scheme.  The  flux-vectors  are  non-linear  functions  of  q,  and  therefore  a  linearization 
procedure  which  retains  the  same  temporal  accuracy  of  eqn.  (17)  has  to  be  introduced 
in  order  to  develop  a  non-iterative  algorithm.  A  linear  equation  with  the  same  temporal 
accuracy  as  eqn.  (17)  can  be  obtained  if  Taylor  series  expansion  is  used  : 

pn+x  ^  pn^  (5«+l  _  5")  +  0[Ae)  (19) 

or,  AF"  =  A"A9"  d- 0( A<2)  (20) 


where 


Similarly, 


where 


A  = 


dq 


{Jacobian  matrix) 


AG"  =  5"A5"  +  0(At=*) 

AVi"  =  {P-R^)^Aq^  +  —{RAq)^  +  0{At^) 

AW2"  =  (g-5„)"Aq"+^(5Ag)"  +  G(At2) 

Orj 


B  = 
P  = 
R  = 


dq 

dVi 

dq 

dVj 

dq( 


5  = 


dW2 

dq 

dW2 

dqr, 


(21) 

(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 
(29) 
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The  constituents  of  the  Jacobian  matrices  are  given  in  Appendix  A. 

A  difficulty  arises  from  the  spatial  cross-derivative  terms  ^(Al/2)”  and  ^(APFi)". 
If  these  terms  were  treated  in  the  same  manner  as  eqn  (19),  considerable  difficulty  will 
be  encountered  in  constructing  an  efficient  spatially  factored  algorithm.  In  order  to 
facilitate  the  approximate  factorization  procedure  it  is  appropriate  and  consistent  to 
treat  the  spatial  cross-derivative  terms  mentioned  above  explicitly  or  eliminate  them 
altogether  in  eqn  (17).  Studies  [2]  on  the  effect  of  the  cross  derivative  terms  in  the 
algorithm  have  indicated  that  there  is  very  little  change  if  the  spatial  cross  derivative 
terms  are  dropped  altogether. 

Substitution  of  eqns.  (20)  to  (24)  into  eqn.  (17),  and  dropping  the  spatial  cross¬ 
derivative  terms  yields 


At 


d  d  9" 


At 


+  Vi  +  v,r+^{-a  +  w,  +  w-,)" 


|A,"  = 
+  0(A«=) 


(30) 


where  the  notations  of  the  form 


d 


and 


^(A  -  p  -h  R,r 


A?"=^[(.4-P  +  JJi)Ajr 

As"  =  A[(b  -Q  +  5,)Asr 


have  been  used  for  convenience. 

Equation  (30)  is  linear,  but  still  would  result  in  a  very  large  matrix  inversion  problem 
when  the  spatial  derivatives  are  replaced  by  finite-difference  formulas.  To  make  the 
problem  simple,  the  left  hand  side  operator  in  eqn.  (30)  is  factored  as  follows  without 
disrupting  the  formal  accuracy  of  the  algorithm. 


1  + 


At 


Ar  R  _  n  +  9  r  _  AA 

dT]^  Q  +  5„)  ^^2 


=  At 


L^(-P  +  -F  F2)"  +  — (-G  +  Wi  +  W2r 


4-  0{At^)  (31) 


The  factored  scheme  (eqn.  (31))  is  implemented  in  the  following  alternating  direction 
fashion 


Am-p-^p  r-  — 

^^(A  P  +  R^)  ^^2 


Aq'=  = 


At 


1^— (-P  -f  Fi  -f  ^2)"  +  -^{-G  -b  +  w^y 


(32) 
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(33) 


^(B-0  +  5,)" 


drj^ 


Ag"  -  Aq^ 


qn+l  ^  ^ 

where  A^^  is  a  dummy  intermediate  vector. 

The  spatial  derivatives  appearing  in  eqns.  (32)  and  (33)  have  to  be  approximated 
by  appropriate  finite  difference  quotients.  Figure  (3)  shows  a  part  of  the  transformed 
uniform  rectangular  grid,  employed  for  solving  the  equations,  where 


6  =  (i-l)A^;  l<i<IL-  r,J  =  (j  -  l)Ar,  l<j<JL  (35) 


where  the  subscripts  i  and  j  refer  to  the  ^  and  rj  direction  respectively.  For  convenience 
A^  and  At/  are  selected  to  be  1.0,  and  hence  have  been  omitted  in  rest  of  the  equations. 
The  transformed  derivatives  x^,  Xr,,  and  ?/^  are  all  computed  using  centered  finite 
differences  at  all  the  interior  points,  and  one-sided,  second-order  accurate  approximations 
along  the  boundaries. 


ri 


j=ji 


i-t-9 

j 

;  1 

Jati 

;  9 

i  ^ 

i- 

2  i- 

1 

i  i- 

hi  i- 

1-2 

^ax.  ^ 


Figure  3:  Transformed  plane. 
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Centered  finite  difference  discretization  is  applied  to  the  spatial  derivatives  in  eqns. 
(32)  and  (33) 


-  W,%.  - 

{I  +  AifeBy  +  Si(-Q  +  5,),^  -  «?5.-,i]“}A,?.  =  Aqtj 


(36) 

(37) 


where  the  finite-difference  operators  used  are 


=  ^i+U 


=  ^ij+h 


(38) 


Beam- Warming  algorithm  is  implemented  in  a  standard  ADI  sequence.  Equation  (36) 
is  solved  first  for  the  vector  Aq‘  for  all  constant  ij  Lines  {2  <  j  <  JL  —  1).  This 
step  involves  the  solution  of  a  block-tridiagonal  linear  system  for  each  constant  rf  line. 
Similarly,  eqn.  (37)  is  solved  next  for  the  vector  Ag”  using  the  available  values  of  the 
intermediate  vector  Aq^  for  each  constant  (  line.  The  vector  q  is  updated  at  time  {n  + 1) 
using  eqn.  (34).  The  updated  q  vector  gives  the  new  flow  variables  at  time  step  (n  4- 1). 


3.2  Numerical  Damping 

A  linear  stability  analysis  of  Beam- Warming  algorithm  indicates  that  the  high  fre¬ 
quency  components  of  the  solution  are  not  damped  when  central  spatial  finite  differences 
are  employed  [55].  Violation  of  stability  condition  produces  an  amplification  of  the  var¬ 
ious  forms  of  error  that  are  present  in  the  numerical  solution.  These  include  truncation 
errors  (due  to  inexact  differentials),  round-off  errors  (due  to  truncated  arithmetics),  and 
errors  due  to  slightly  inconsistent  boundary  conditions.  The  solution  becomes  stable 
when  some  artificial  dissipation  is  added  to  the  algorithm.  The  artificial  dissipation 
term  also  provides  smoothing  for  the  capturing  of  embedded  shocks.  The  present  study 
employs  fourth  order  explicit  and  second  order  implicit  dissipation  terms  [32].  The 
fourth  order  explicit  dissipation  term  removes  the  small  wavelength  oscillations  (2Aa: 
oscillation)  in  the  flowfield.  The  second  order  implicit  dissipation  term  actually  stabi¬ 
lizes  the  system  by  increasing  the  diagonal  dominance  of  the  block  tridiagonal  matrix 
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formed  in  the  algorithm.  The  fourth-order  explicit  artificial  dissipation  term  is  appended 
to  the  right  hand  side  of  eqn.  (36).  The  implicit  dissipation  term  is  added  to  the  left 
hand  side  of  eqn.  (36).  The  addition  of  the  artificial  dissipation  terms  results  in  the 
following  algorithm  : 


{i+Y  "}  +  OKi 

-Ai[,.;(F  -  +  ft,(G  -  W,),.,.  -  (39) 


{/  +  AtUijBij  +  S,i-Q  +  S,)y  -  =  Ag; 


(40) 


where  Wj  and  We  are  the  coefficients  of  implicit  and  explicit  artificial  dissipation  terms, 
respectively.  Addition  of  the  second-order  implicit  damping  term  extends  the  linear 
stability  bound  of  the  fourth-order  explicit  damping  term. 


3.3  Geometric  Conservation  Law 


The  numerical  simulation  of  unsteady  flow  past  a  moving  airfoil  involves  the  move¬ 
ment  of  the  computational  grid.  Maintenance  of  global  conservation  of  mass,  momentum, 
and  energy  is  an  important  part  of  numerical  simulation.  Use  of  boundary-conforming 
coordinate  transformations  and  the  subsequent  application  of  finite-difference  formulas 
to  the  moving  boundaries  or  grids  lead  to  grid-movement  related  errors  in  the  solution. 
To  remove  the  grid-movement  related  errors,  the  Geometric  Conservation  Law  (GCL)  is 
implemented  in  the  algorithm  [44].  Details  of  Geometric  Conservation  Law  is  presented 
in  Appendix  B.  The  following  term  is  added  to  the  right  hand  side  of  the  algorithm  to 
satisfy  the  global  conservation  (GCL). 


GCL  term  =  At 


(41) 


where 


6  =  (-®ty*?  +  ytXr,)J ;  vt  =  {xty^  -  ytXi)J 


(42) 
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3.4  Computational  Grid 


Nearly  orthogonal  boundary-fitted  grids  have  been  employed  for  the  present  study. 
For  numerical  simulation  of  flow  past  an  airfoil,  three  types  of  boundary-fitted  grids 
can  be  used  :  0-grid,  C-grid  and  H-grid.  In  the  present  study,  C-type  grids  have  been 
utilized  for  the  simulation.  An  advantage  of  the  C-grid  is  improved  grid  resolution 
near  the  trailing  edge.  A  schematic  diagram  of  a  C-grid  has  been  shown  in  Fig.  (4). 
The  grids  employed  in  the  present  study  have  been  generated  by  the  hyperbolic  grid 
generation  code  developed  by  Kinsey  and  Barth  [19].  The  grid  generation  code  needs  the 
distribution  of  points  on  the  airfoil  surface  and  wake  as  input  and  gives  the  coordinates 
of  the  nodes  of  the  C-grid  as  output.  Figure  (5)  shows  a  part  of  the  C-grid  near  the 
airfoil. 


Figure  4:  C-grid  topology  and  computational  boundaries. 
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Figure  5:  C-grid  for  NACA-0012  airfoil. 

3.5  Boundary  Conditions 

Suitable  boundary  conditions  have  to  be  specified  to  completely  define  the  problem. 
The  boundary  conditions  are  based  on  the  particular  problem  being  considered.  In  the 
present  study,  a  C-grid  has  been  employed.  A  C-grid  has  numerical  boundaries  at  the 
airfoil  surface,  the  outer  boundary  (inflow  and  outflow),  and  C-grid  cut  in  the  wake 
region.  The  boundary  condition  at  all  boundaries  except  at  the  C-grid  cut  in  wake 
region  are  imposed  explicitly. 

On  the  airfoil  surface,  the  following  adiabatic,  no-slip  condition  is  applied  : 

tl  =  UB  (43) 


The  normal  pressure  gradient  incorporates  the  effect  of  the  airfoil  motion. 


=  —pobB  • 
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where  Ub  and  as  are  the  velocity  and  acceleration  of  the  airfoil  surface,  respectively. 
For  the  rotation  of  the  airfoil  about  a  fixed  axis  with  a  pitch  rate  Ub  and  ob  can  be 
defined  as  : 


^  X  (tb 


ro) 


dB  =  ^  X  [rB 


fo)  -  ^^{tb  -  ro) 


(46) 

(47) 


To  is  the  position  vector  of  the  point  of  rotation  and  tb  is  the  position  vector  of  the 
point  on  the  airfoil  where  the  velocity  and  the  acceleration  has  to  be  determined. 

At  the  infiow  and  outflow  boundaries,  a  one- dimensional  Method  of  Characteristics 
boundary  condition  [45]  is  applied.  The  Method  of  Characteristics  is  described  in  de¬ 
tail  in  Appendix  C.  A  local  coordinate  system  (x,  y)  is  constructed  orthogonal  to  the 
outer  boundary  with  x  normal  to  and  directed  outward  from  the  boundary  and  y  di¬ 
rected  along  the  boundary.  Assuming  that  the  derivatives  along  the  boundary  can  be 
neglected  (one-dimensional),  the  following  characteristic  equations  and  corresponding 
characteristic  variables  have  been  found  : 

2d  dx 

Characteristic  equation  1  (C^)  :  u-\ - -  =  constant  along  —  =  u  a 

7  —  1  at 


Characteristic  equation  2  (C^) 


_  2a 

u - =  constant 

7-1 


dx  _ 

along  =  u  —  a 
dt 


(49) 


Characteristic  equation  3 


r) 


s  =  constant 


along 


dx 


(50) 


Characteristic  equation  4  (C**)  :  v  =  constant  along  —  =  u  (51) 

dt 

where  s  is  the  entropy,  a  is  the  speed  of  sound,  and  u  and  v  are  velocities  in  x  and  y 
directions  respectively. 

All  of  the  Riemann  variables  corresponding  to  the  characteristics  moving  out  of  the 
computational  domain  have  to  be  interpolated  from  the  values  available  inside  the  do¬ 
main.  All  other  Riemann  variables  are  set  as  boundary  conditions.  The  outer  boundary 
is  defined  at  a  sufficiently  large  distance  from  the  airfoil  to  insure  accuracy  of  this  bound¬ 
ary  condition  (see  section  4.2). 
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(i)  Transformed  domain 


(ii)  Wake  region  of  computational  domain 


Figure  6:  Computational  and  corresponding  transformed  domain  of  the  region  near  the 
trailing  edge  showing  the  usage  of  periodic  condition  at  the  C-grid  cut. 
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Figure  7:  Transformed  domain  during  the  sweep. 

At  the  C-grid  cut  (wake  region),  the  flow  variables  are  solved  implicitly  by  imposing 
the  periodicity  condition.  Figure  (6)  shows  the  region  near  the  trailing  edge  and  the 
transformed  domain.  The  transformed  domain  is  obtained  by  unwrapping  the  compu¬ 
tational  grid  from  around  the  airfoil.  In  the  transformed  domain  the  value  of  the  flow 
variables  at  point  ‘a’  are  the  same  as  at  point  ‘d’.  Similarly,  values  at  ‘m’  and  ‘n’  are 
equal,  and  at  ‘6’  and  ‘c’  (trailing  edge)  are  equal  and  so  on.  This  periodicity  condition 
is  used  to  solve  for  the  flow  variables  implicitly  at  the  C-grid  cut.  The  Beam- Warming 
algorithm  is  solved  utilizing  the  ADI  scheme.  The  transformed  domain  looks  like  Fig.  (6 
(i))  during  the  ^-sweep.  But  during  the  7/-sweep,  the  transformed  domain  looks  like 
fig.  (7).  This  way  the  variables  at  the  C-grid  cut  can  be  solved  implicitly. 
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4  GOVERNING  EQUATIONS 
(UNSTRUCTURED  GRID) 


The  governing  equations  are  the  2-D  compressible  laminar  Navier- Stokes  equa¬ 
tions.  For  an  arbitrary  deformable  control  volume  of  volume  (per  unit  depth)  and 
surface  dV ^  the  nondimensional  equations  are 

/  Qdxdy  -I-  /  {Fdy  -  Gdx)  =  0  (52) 

dt  Jv  JdV 

where  Q  is  the  vector  of  dependent  variables 

Q  =  {p,pu,pv,pef  (53) 

where  p  is  the  density,  u  and  v  are  the  velocity  components  in  the  x—  and  y— directions, 
and  e  is  the  total  energy  per  unit  mass, 

"  =  .,(.^-1)+^""+’’")  (54) 

where  T  is  the  static  temperature  and  7  =  c^/c„  is  the  ratio  of  specific  heats.  The  flow 
variables  are  nondimensionalized  using  the  dimensional  reference  quantities  L  (length), 
Poo  (density),  Uco  (speed  of  sound)  and  Too  (static  temperature).  The  superscript  “t” 
denotes  the  vector  transpose.  The  flux  vectors  are 

pU 

puU  +  p  —  Txx 
pvU  -  T^y 
peU  +  jm  +  Px 
pV 

puV  ^xy 
PVV  +  P-Tyy 

peV  +  pv  +  j3y 
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where 


U  =  u  —  Ug,  V  —  V  —  Vg 


(57) 


where  Ug  and  Vg  are  the  x—  and  ?/— components  of  the  velocity  of  the  surface  of  the 
control  volume.  The  components  of  the  viscous  stress  tensor  are 


'^XX 


Txy 


‘yy 


Moo  {^du  2dv\ 
Reoo^\3dx  Sdy) 
Moo  f  9u 
Reoo^  \dy  dxj 

Moo 

R^oo^  \3  dy  3  dx / 


(58) 


where  Moo  is  the  reference  Mach  number,  Reoo  is  the  reference  Reynolds  number,  and 
/j,  is  the  molecular  viscosity  normalized  by  a  suitable  reference  value  fioo-  The  static 
pressure  p  is  normalized  by  and  satisfies  the  Ideal  Gas  Equation 


PT 


Furthermore, 


(59) 


—  Qx  '^xx^  '^xy^ 

Py  ~  Qy  ~  '^xy'U'  —  TyyV  (60) 


where  the  heat  flux  is 


Moo  1  dT 
^  Rcoo  Pr  (7— 1)  dx 
Moo  1  dT 

^  Reoo  Pri-y-l)  dy 


(61) 


where  Pr  is  the  Prandtl  number. 

It  is  convenient  to  decompose  the  flux  vectors  into  their  inviscid  and  viscous  contri¬ 
butions 

pU 

puU  d-  p 

Finv  = 

pvU 

pell  -f  pu 
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0 


and 


F  ■  =  < 

Via  —  ' 


— r. 


xy 


13. 


Gir..  =  < 


Gi,i»  —  < 


pV 
puV 
pvV  +  p 
peV  +  pv 

0 


—T, 


xy 


—T, 


yy 


(62) 


(63) 
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NUMERICAL  ALGORITHM 
(UNSTRUCTURED  GRID) 


5.1  Implicit  Algorithm 

An  unstructured  grid  of  triangles  is  employed  (Fig.  8).  A  cell-centered  storage  archi¬ 
tecture  is  assumed  where  i  denotes  the  cell  index.  The  cell-averaged  value  of  the  vector 
of  dependent  variables  Q  is 

Qi  =  ^^Jv  Qdxdy  (64) 

where  Vi  is  the  volume  (per  unit  depth)  of  the  cell.  The  governing  equations  (52)  are 
therefore 

f^{QiVi)  +  Ci  =  0  (65) 

where  Ci  is  the  net  flux  across  the  cell  faces 


Ci=  [  (Fdy  -  Gdx) 

JdVi 


(66) 


It  is  furthermore  assumed  that  Vi  is  constant,  although  this  restriction  can  be  easily 
relaxed. 

A  family  of  implicit  algorithms  for  eqn.  (65)  is 


/  + 


l+aVidQij  1  +  aV, 


ie* 

1  At _ _  a 


c'r + 


n-l 


1  -f  a  Vi  ’  1  -f  a 


where  At  is  the  time  increment,  and 


(67) 


AQ?  =  Q"«  -  Q: 


(68) 


and  denotes  the  summation  over  all  cells  constituting  the  numerical  domain  of 

dependence  (“star”)  in  the  application  of  eqn.  (65)  at  cell  i  with  the  exception  of  the 
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Figure  8:  Unstructured  grid  of  triangles 


Table  1:  Family  of  Implicit  Algorithms 


Type 

Accuracy 

a 

Euler 

0(At) 

0 

1 

Trapezoidal 

0(Ai2) 

0 

1 

2 

3-point 

C>(Ai2) 

1 

2 

1 

cell  i.  Note  that  the  Einstein  summation  convention  is  not  employed  in  eqn.  (67).  The 
family  of  algorithms  and  the  associated  temporal  accuracy  is  indicated  in  Table  1. 

The  increments  A(5i  for  the  N  cells  can  be  concatenated 

AQ  =  (AQi,A(^2,...,Ag;,,)^  (69) 

based  on  an  assumed  ordering  of  the  cells.  Application  of  eqn.  (67)  to  each  cell,  together 
with  suitable  boundary  conditions,  yields  the  linear  system  of  equations, 

A^,Q  =  n  (70) 

where  A  is  the  generalized  Jacobian  matrix,  and  TZ  is  the  residual.  The  Jacobian  A  is 
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sparse  and  banded. 


5.2  Inviscid  Fluxes 


The  contribution  to  Ci  from  the  inviscid  fluxes  is  obtained  from  flux  difference  split 


method  of  Roe  [34] 


fc=3 


-  dx)  =  Y.  T-'^HAs 
fc=l 


where  the  rotation  matrix  is 


(71) 


Ax  0  0  0 


T-'^As  =  \ 


0 

0 


Ay  Ax  0 
—Ax  Ay  0 


(72) 


0  0  0  As 


where  Ax  and  Ay  are  the  x—  and  y—  projections  of  side  k  of  cell  i  and  As^  —  Ax^  +  Ay^. 
The  flux  vector  is 


pU 

pUu  +  p 
pvU 

peU  +  pu 


(73) 


where 


u 


V 

U 


Ay  Ax 

u- - V— 

As  As 


Ax 

u— - 1-  V 

As 


Ay 

As 


U^-V— 

As  As 


(74) 


The  flux  vector  if  on  a  face  fc  of  a  triangle  is  approximated  by 


if  =  I  (ffj  +  ii,  +  S\k\S-^AR)  (75) 

The  terms  Hi  and  Hr  denote  the  “left”  and  “right”  reconstructed  values  of  the  flux 
vector  corresponding  to  the  inside  surface  {i.e.,  inside  the  triangle)  and  outside  surface 
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of  face  k.  These  values  are  obtained  by  a  spatially  second-order  accurate  reconstruction 
of  Q  using  the  cell-centered  values  [57] 


Q  =  Qi-^  •  r  (76) 

where  Qi  denotes  the  value  of  Q  at  the  cell  centroid,  f  is  the  vector  from  the  cell  centroid 
to  the  midpoint  on  the  face,  and  $  is  a  limiter  function  [57]  which  is  determined  by  the 
requirement  that  the  reconstructed  value  of  Q  at  the  cell  face  is  within  the  range  of 
values  of  Q  for  those  cells  employed  in  the  reconstruction.  The  gradient  VQ  at  the  cell 
centroid  is  obtained  from  Green’s  Theorem  using  the  dual  triangle  abc  whose  vertices 
are  the  centroids  of  the  cells  adjacent  to  cell  i  (Fig.  9).  Thus, 

Q  =  Qi  ^  (AQa  +  AQb  -|-  AQc)  (77) 


Figure  9:  Triangle  used  to  construct  VQ 
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The  additional  term  in  eqn.  (75)  represents  the  contribution  to  the  flux  vector  from 
the  waves  originating  within  the  adjacent  cells  and  may  be  written, 

S|A|S-'AB  =  2a,|A^|e,-  (79) 

J=1 

where  |Aj|  are  the  absolute  values  of  the  eigenvalues 

Ai  =  ^ 

k  =  u 

A3  =7  U  a 

A4  -  fi-a  (80) 

The  Roe-averaged  variables  are 

^  "b 

~  _  y/plUl  +  y/^Ur 

y/^  +  y/^ 
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y/Fl  +  y/^ 

^  _  y/^l  Hi  -|-  y'pjT-Rr 
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=  (7  -  1)  (H  -  q) 

9  =  +  (81) 

The  eigenvectors  are 
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The  coefficients  are 
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where,  for  example,  Apu  =  pui  —  pur. 

An  entropy  cut-olf  is  employed  to  eliminate  unphysical  expansion  shocks  which  can 
occur  when  one  or  more  of  the  eigenvalues  are  zero.  In  eqn.  (79),  |Aj|  is  replaced  by 
II  Xj  II  where 


|A.- 

(|A,-P  +  5=^)  /28 


if  |A,-|  >  6 
if  |A,-|  <  6 


(84) 


where  5  is  a  small  quantity  (typically,  S  <  0.1). 

The  individual  elements  of  the  additional  flux  term  ^'lAl.S'^Ai?  are  given  in  Ap¬ 
pendix  D. 


5.3  Viscous  Fluxes 

The  contribution  to  Ci  from  the  viscous  fluxes  and  heat  transfer  on  face  k  is  obtained 
from  application  of  Gauss’  Theorem  [22]  to  the  quadrilateral  defined  by  the  cell  centroids 
of  the  cells  adjacent  to  face  k  and  the  two  nodes  defining  the  endpoints  (Fig.  10).  For 
any  function  f{x,y), 
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where  V  and  dV  are  the  volume  and  surface  of  the  quadrilateral  abed,  respectively,  and 
Tlx  is  the  component  of  the  outwards  normal  in  the  x— direction.  A  similar  equation  is 
obtained  for  df  jdy.  The  molecular  viscosity  is  evaluated  at  the  midpoint  of  surface  k 
using  the  formula 

/i  =  K/ife  +  fid)  (86) 

where  fii  and  fij  represent  the  molecular  viscosity  evaluated  at  nodes  b  and  d,  respec¬ 
tively. 

The  values  of  Q  at  the  nodes  are  needed  only  for  the  viscous  terms  and  are  obtained 
by  second-order  interpolation  of  Q  from  those  cells  sharing  the  node  [33] 


Qi  =  £  Yj 

cells  cells 


(87) 


where  Qj  denotes  the  Q  at  node  j,  Qi  denotes  Q  at  the  centroid  of  cell  i  which  shares 
the  node  j,  the  sum  is  over  all  cells  sharing  the  node  {xj,yj),  and  Wi  are  dimensionless 
weights 
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5.4  Jacobian  of  Inviscid  Terms 

The  exact  expression  for  the  individual  Jacobians  of  the  inviscid  terms  in  eqn.  (67) 
is 
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(90) 
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Figure  10:  Quadrilateral  employed  for  determination  of  viscous  fluxes  and  heat  transfer 
where 


A,  = 

A.  = 

Bi  = 

ag, 

Br  = 

FiC 

C  = 

|r-^5|A|5-^rA5 

(91) 


Complete  details  of  the  exact  inviscid  Jacobian  are  provided  in  Appendix  D. 


5.5  Jacobian  of  Viscous  Terms 


The  exact  expression  for  the  individual  Jacobians  of  the  viscous  terms  in  eqn.  (67)  is 
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Complete  details  of  the  exact  viscous  Jacobian  are  provided  in  Appendix  D. 

5.6  Boundary  Conditions 

Boundary  conditions  are  incorporated  implicitly.  At  a  supersonic  inflow  boundary, 
Q  is  specifled  fully.  At  a  solid  boundary,  the  velocity  is  set  equal  to  zero,  the  surface 
temperature  or  heat  flux  {e.g.,  adiabatic)  is  specified,  and  the  normal  derivative  of  the 
static  pressure  is  set  to  zero.  On  a  symmetry  boundary,  the  normal  component  of  the 
velocity  is  set  to  zero,  and  the  normal  gradients  of  the  remaining  flow  variables  are  set 
to  zero.  At  a  supersonic  outflow  boundary,  the  flow  variables  are  extrapolated  linearly 
from  the  interior. 

5.7  Solution  of  Linear  System 

The  linear  system  eqn.  (70)  is  solved  at  each  time  step  using  the  BiCGSTAB  algorithm 
of  Van  der  Vorst  [46].  The  algorithm  is  a  conjugate  gradient  type  iterative  method. 
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6  CODE  VALIDATION  (STRUCTURED  GRID) 


A  two-dimensional,  laminar,  compressible  Navier-Stokes  solver  has  been  devel¬ 
oped  employing  the  Beam- Warming  algorithm  described  in  section  3.  A  variety  of  two- 
dimensional  test  computations  were  performed  to  establish  the  accuracy  of  the  solver. 
The  solver  was  validated  through  application  to  the  following  computations  : 

(а)  Flat  plate  boundary  layer. 

(б)  Steady  state  laminar  flow  past  a  NACA-0012  airfoil. 

(c)  Unsteady  laminar  flow  past  a  pitching  NACA-0015  airfoil. 

6.1  Flat  Plate  Boundary  Layer 

The  supersonic  laminar  boundary  layer  on  a  flat  plate  provides  a  test  of  the  accuracy 
of  the  solver  for  viscous  flows  by  comparison  with  an  exact  solution  of  the  boundary 
layer  equations  under  the  conditions  Pr  =  1,  and  fi  =  T.  The  solver  was  employed  to 
perform  a  boundary  layer  computation  over  a  flat  plate.  The  freestream  Mach  number 
Moo  was  2.0  and  Reynolds  number  based  on  the  length  of  the  plate  {Rcl)  was  selected 
to  be  10^. 

The  computational  configuration  is  shown  in  Fig.  (11).  The  inflow  is  a  uniform 
steady  flow  and  extrapolation  was  used  to  find  the  flow  variables  at  the  outflow.  The 
Method  of  Characteristics  was  used  at  the  top  boundary.  At  the  lower  boundary,  in 
front  of  the  plate  leading  edge,  a  symmetry  boundary  condition  was  used  for  the  flow 
variables  in  the  7/-direction  (i.e.,  )  =  0).  At  the  plate  boundary,  a  no-slip,  adiabatic 

boundary  condition  was  used  for  the  velocity  and  temperature,  respectively.  A  linear 
viscosity  law  was  employed  with  Pr  =  1.0. 

JL  =  1- 

/^OO  I'oo 

where  the  subscript  oo  denotes  the  freestream  values. 
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top  boundary 


Figure  11:  Schematic  of  the  computational  domain  for  boundary  layer  over  a  flat  plate. 


The  grid  employed  =  60  and  Nr,  —  75  points  in  the  streamwise  (^)  and  normal 
{rj)  directions,  respectively.  The  number  of  grid  points  within  the  boundary  layer  was 
Nbl  =  54  at  the  end  of  the  flat  plate. 

The  computed  velocity  and  static  temperature  profiles  are  compared  with  the  avail¬ 
able  analytical  profiles  [56]  in  Figures  (12)  and  (13)  at  x  =  0.8T.  The  ordinate  is  the 
transformed  distance  in  rj  direction: 


_p_^ 
poo  X 


(95) 


The  velocity  is  normalized  by  the  freestream  velocity,  Uoo,  and  the  static  temperature 
is  normalized  by  the  freestream  static  temperature.  Too-  The  agreement  between  the 
computation  and  the  theory  is  excellent.  The  computed  velocity  and  static  temperature 
profiles  agreed  with  the  analytical  profiles  to  within  1.0%. 

A  computation  was  also  performed  with  the  grid  at  an  angle  to  the  horizontal. 
The  results  showed  similar  agreement  with  the  theory  as  in  the  previous  case.  This 
computation  ascertained  the  accuracy  of  the  grid  transformation  from  (x,  y)  coordinates 
to  coordinates. 
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6.2  Steady  Flow  Past  a  Stationary  NACA-0012  Airfoil 


The  viscous  subsonic  flow  past  an  airfoil  at  a  fixed  angle  provides  a  realistic  test  of 
the  accuracy  of  the  numerical  algorithm  by  comparison  with  previous  numerical  com¬ 
putations. 

A  series  of  computations  were  performed  for  a  stationary  NACA-0012  airfoil  at  an 
angle  of  attack,  a  =  0°,  Reynolds  number  based  on  the  chord  Rcc  =  5  x  10^,  and 
freestream  Mach  number  Moo  =  0.2.  Sutherland’s  Law  was  employed,  and  Pr  =  0.73.  A 
C-grid  was  employed  with  =  303  and  Nr,  =  101,  where  ^  is  the  curvilinear  coordinate 
in  the  direction  of  the  “(7”,  and  1]  is  the  coordinate  approximately  orthogonal  to  rj. 
The  average  normal  distance  of  the  first  row  of  mesh  points  adjacent  to  the  airfoil  was 
Auauer  =  5.0  X  10“^c,  and  25  grid  points  were  contained  within  the  boundary  layer  at 
mid-chord.  The  grid  spacing  tangential  to  the  airfoil  surface  varied  from  As  =  4.3  x  10“'* c 
to  1.99  X  10“^c,  with  the  finest  grid  employed  near  the  leading  and  trailing  edges. 

The  results  were  compared  with  previous  computations  by  Reran  [3]  and  the  results 
utilizing  the  present  unstructured  grid  algorithm  [20]  at  the  same  flow  field  conditions, 
and  Mehta  [26]  at  the  same  Reynolds  number  for  incompressible  flow.  Table  2  presents 
the  computed  drag  coefficient  of  the  present  study  and  the  three  other  computations. 
Figure  14  shows  the  comparison  of  the  computed  surface  pressure  coefficient  with  Re¬ 
run’s  results.  The  computed  drag  coefficient  and  surface  pressure  coefficient  show  close 
agreement  with  the  other  computations. 

Table  2:  Computed  Drag  Coefficient  (Stationary  NACA-0012  Airfoil),  Rcc  =  5  x  10^, 
Moo  =  0.2,  a  =  0°. 


Reference 

drag  coefficient  [cj) 

Present  study  (Structured  grid) 

0.0544 

Reran 

0.0530 

Present  study  (Unstructured  grid) 

0.0547 

Mehta 

0.0534 
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Figure  14:  Pressure  coefficient  on  the  surface  of  a  stationary  NACA-0012  airfoil  at  a  =  0° 
{Re,  =  5  X  10^  Moc  =  0.2). 

The  Beam- Warming  algorithm  has  artificial  damping  terms.  It  is  important  to  study 
the  effect  of  the  damping  coefficients  on  the  computations.  A  number  of  computations 
were  performed  to  study  the  sensitivity  of  the  computed  flowfield  to  the  magnitude  of  the 
explicit  damping  coefficient,  u),.  The  steady  solutions  do  not  depend  on  the  magnitude  of 
the  implicit  dissipation  terms  because  at  steady  state  the  right  hand  side  of  equations  39 
and  40  are  zero  and  the  implicit  dissipation  term  does  not  have  any  effect  on  the  final 
solution.  The  computed  flowfield  was  found  to  be  insensitive  to  the  magnitude  of  the 
explicit  damping  coefficient  u},  provided  u;,  <  0.01.  Figure  (15)  shows  the  values  of 
the  drag  coefficient  and  leading  edge  pressure  coefficient  for  different  values  of  explicit 
damping  coefficient. 

The  effect  of  the  location  of  the  outer  boundary  distance  from  the  airfoil  on  the  com¬ 
puted  flowfield  was  also  studied.  Figure  (16)  shows  that  the  computed  drag  coefficient 
and  leading  edge  pressure  coefficient  are  insensitive  to  the  location  of  the  outer  boundary 
provided  the  distance  from  the  airfoil  exceeds  18c. 
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Figure  15:  Effect  of  explicit  damping  coefficient,  Wg,  on  the  drag  coefficient  and  leading 
edge  pressure  coefficient  for  a  stationary  NACA-0012  airfoil  at  a  =  0°  {Rtc  =  5  x  10®, 

Moo  =  0.2). 


6.3  Unsteady  Flow  Past  a  Pitching  NACA-0015  Airfoil 

The  subsonic  viscous  unsteady  flow  past  a  pitching  airfoil  provides  a  test  of  the 
temporal  and  spatial  accuracy  of  the  numerical  algorithm  through  comparison  with 
previous  numerical  simulations. 

A  set  of  computations  were  performed  to  study  the  flow  past  an  NACA-0015  airfoil 
pitching  about  its  quarter  chord  axis  at  Re^  =  10“*,  M^o  =  0.2,  and  non-dimensional 
pitch  rate  =  0.2.  The  pitching  motion  was  defined  by  : 

n+{t)  =  (^1  -  (96) 

where  to  is  the  time  at  which  the  non-dimensional  pitch  rate  has  reached  99%  of  its 
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outer  boundary  distance  in  chord  lengths 


Figure  16:  Effect  of  the  location  of  the  outer  boundary  distance  from  the  airfoil  on  the 
drag  coefficient  and  leading  edge  pressure  coefficient  for  a  stationary  NACA-0012  airfoil 
at  a  =  0°  (iZcc  =  5  X  10^  Moo  =  0.2). 

asymptotic  value,  The  non-dimensional  pitch  rate  is  defined  as  : 

(97) 

y  oo 

where  fl!  is  the  pitching  rate  in  radfsec,  c  is  the  chord  length  and  Uoo  is  the  freestream 
velocity.  The  functional  form  of  the  non-dimensional  pitch  rate  provides  a  smooth 
acceleration  of  the  airfoil  to  its  asymptotic  value  during  an  effective  time  interval  to  = 
O.bc/Uoo  which  corresponds  to  4.57°  of  pitch.  The  functional  form  of  the  pitching  motion 
also  removes  the  infinite  acceleration  faced  at  the  start  of  the  pitching  if  a  constant 
value  of  is  applied  throughout  the  motion.  The  flow  condition  and  expression  for 
the  pitching  motion  correspond  to  the  previous  computation  by  Visbal  [49]. 

A  C-grid  was  employed  with  =  332  and  =  85.  The  average  normal  distance 
of  the  first  row  of  mesh  points  adjacent  to  the  airfoil  was  Auaver  =  1-0  x  10“^c,  and  33 
grid  points  were  contained  within  the  boundary  layer  at  mid-chord.  The  grid  spacing 
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tangential  to  the  airfoil  surface  varied  from  As  =  4.59  x  10~®c  to  1.47  x  10~^c,  with  the 
finest  grid  employed  near  the  leading  and  trailing  edges. 


Figure  17:  Lift,  drag  and  moment  coefficients  for  a  pitching  NACA-0015  airfoil  (iJcc  = 

10^  Moo  =  0.2,  =  0.2). 

The  computed  flowfields  were  found  to  be  in  close  agreement  with  the  previous 
computation  by  Visbal  [49],  who  used  an  0-grid.  Figure  (17)  shows  the  comparison  of 
the  computed  lift  coefficient  (7j,  drag  coefficient  Cd,  and  moment  coefficient  Cm  with  the 
previous  results  of  Visbal  for  0  <  a  <  25°,  where  a  is  the  angle  of  attack  of  the  airfoil. 
Detailed  examination  of  the  flow  variables  showed  similar  close  agreement. 
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7  CODE  VALIDATION 
(UNSTRUCTURED  GRID) 


7.1  Riemann  Shock  Tube 

The  Riemann  Shock  Tube  Problem  provides  a  test  of  the  temporal  and  spatial  ac¬ 
curacy  of  the  numerical  algorithm  for  inviscid  flows  through  comparison  with  an  exact 
solution  [21].  A  total  of  four  rows  of  200  triangular  cells  each  were  employed  with  cell 
spacing  Ax  =  Ay  =  1.0.  The  pressure  ratio  p^/pi  =  10,  where  p^  and  pi  are  the  static 
pressure  in  the  driver  and  driven  sections,  respectively.  The  static  temperature  ratio 
T4/T1  =  1.0.  The  diaphragm  is  located  in  the  middle  of  the  computational  domain.  The 
computations  employed  trapezoidal  time  differencing  with  a  Courant  number  of  2.  The 
entropy  cut-oflf  8  =  0.1. 

Results  for  the  static  pressure,  static  temperature  and  velocity  are  shown  in  Figs.  18 
through  20  at  t  =  18.4.  The  computed  profiles  are  in  close  agreement  with  the  exact 
result  [21]. 

7.2  Flat  Plate  Boundary  Layer 

The  supersonic  laminar  boundary  layer  on  a  flat  plate  provides  a  test  of  the  spatial 
accuracy  of  the  numerical  algorithm  for  viscous  flows  through  comparison  with  an  exact 
solution  of  the  boundary  layer  equations  under  the  conditions  Pr  =  1  and  p  =  T  [56]. 
A  total  of  7802  triangular  cells  were  utilized  comprising  47  rows  of  166  cells  each.  The 
streamwise  cell  length  Ax  =  1.28  x  10“^.  The  transverse  ceU  length  Ay  varied  from 
1.0065  X  10“®  adjacent  to  the  flat  plate  to  2.07  x  10“^  at  the  upper  boundary  of  the 
computational  domain.  The  Reynolds  number  Rcoo  based  on  plate  length  was  3.95  x  10^. 

The  computed  and  theoretical  velocity  and  static  temperature  profiles  are  displayed 
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Temperatur* 


Figure  18:  Comparison  of  static  pressure  in  a  shock  tube 


Figure  19:  Comparison  of  static  temperature  in  a  shock  tube 
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1.25 


Figure  20:  Comparison  of  velocity  in  a  shock  tube 
in  Figs.  21  and  22,  respectively.  The  ordinate  is  the  transformed  dimensionless  distance 

’l  =  Y^,Y  =  J%dy  (98) 

The  velocity  is  normalized  by  the  freestream  speed  of  sound  (and  hence  u^o  =  Moo), 
and  the  static  temperature  is  normalized  by  the  freestream  static  temperature.  The 
agreement  between  the  computation  and  theory  is  excellent.  The  velocity  profile  agrees 
with  the  theoretical  result  to  within  1.5%.  Similar  close  agreement  is  observed  with  the 
static  temperature  profile.  In  particular,  the  computed  wall  temperature  agrees  with 
the  theoretical  value  {T^  =  1.80)  within  0.5%. 

7.3  NACA-0012  Airfoil 

The  subsonic  viscous  fiow  past  an  NACA-0012  airfoil  at  a  =  0  deg  was  computed  and 
compared  with  previous  numerical  solutions.  A  total  of  25,534  cells  were  utilized.  The 
ceU  nodes  were  generated  using  the  GRAPE  program  of  Sorenson  [42]  and  the  Delaunay 
triangulation  obtained  using  the  UNSTRUCT  program  of  Merriam  [27].  The  nodes  adja¬ 
cent  to  the  airfoil  surface  were  typically  at  a  normal  distance  of  2.0  x  10“^c  where  c  is 
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Figure  21:  Comparison  of  velocity  in  the  boundary  layer 


Figure  22:  Comparison  of  static  temperature  in  the  boundary  layer 
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the  airfoil  chord.  The  freestream  Mach  number  Moo  —  0.2  and  the  Reynolds  Number 
Rcc  =  5  X  10®  based  on  the  chord. 

The  computed  drag  coefficient  is  compared  with  the  structured  grid  results,  results 
due  to  Reran  [3],  and  the  results  of  Mehta  [26]  (Table  2).  The  flow  conditions  of  all 
computations  are  identical,  except  the  Mach  number  was  zero  for  the  computations  by 
Mehta.  The  computed  values  of  cj  are  in  close  agreement. 


X/C 

Figure  23:  Comparison  of  surface  pressure  coefficient  on  a  NACA-0012  airfoil 

The  computed  surface  pressure  is  shown  in  Fig.  23  for  the  present  algorithm,  together 
with  the  results  of  structured  grid  computation  and  results  of  Reran  [3]  which  both 
employed  the  Ream- Warming  algorithm.  The  agreement  is  excellent.  The  computed 
leading  edge  pressure  coefficient  Cp  agrees  with  the  exact  (inviscid)  value  (cp  =  1.01) 
Avithin  2%. 
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8  DEFINITION  OF  PROBLEM 


The  research  problem  is  described  in  the  present  section.  A  brief  description  of 
the  critical  point  theory  follows  the  definition  of  the  problem.  The  critical  point  theory 
has  been  applied  to  the  computation  of  unsteady  separation  on  a  pitching  airfoil  to 
better  understand  the  formation  and  development  of  the  structures  in  the  fiowfield. 

8.1  Problem  Description 

The  focus  of  the  research  is  understanding  of  the  effects  of  compressibility,  pitch 
rate  and  Reynolds  number  on  the  incipient  boundary  layer  separation  for  viscous  flow 
past  a  NACA-0012  airfoil  pitching  about  the  quarter  chord  axis.  The  flow  and  C-grid 
configurations  are  shown  in  Figure  (24).  The  pitching  motion  is  defined  by  equation 
(96).  The  pitching  motion  was  initiated  after  the  fiowfield  had  been  fully  established  at 
a:  =  0°.  The  initial  fiowfield  displayed  a  slight  unsteadiness  due  to  the  periodic  vortex 
shedding  at  the  trailing  edge.  All  computations  were  initiated  when  the  lift  coefficient 
was  zero  and  the  results  presented  herein  are  insensitive  to  the  relative  phase  of  the 
initiation  of  the  pitching  motion. 

An  extensive  grid  refinement  study  was  conducted  for  each  of  the  cases  studied  uti¬ 
lizing  the  structured  grid  algorithm  to  establish  the  accuracy  of  the  solutions.  Complete 
details  of  the  grids  employed  for  the  different  cases  are  given  in  Table  3.  A  total  of  7 
cases  were  studied  (figure  2).  The  accuracy  of  the  computations  were  assessed  for  each 
case  based  on  the  comparison  of  the  results  for  two  grids  with  successive  refinement.  A 
non-dimensional  time  step  {df^  =  tUoolc)  of  1.0  x  10“®  was  employed  for  all  the  seven 
cases  studied  utilizing  the  structured  grid  algorithm.  Computations  conducted  with  a 
time  step  of  0.5  x  10“^  demonstrated  that  the  computed  results  were  insensitive  to  the 
selected  time  step. 

One  computation  was  performed  for  Case  1  using  the  unstructured  grid  algorithm 
for  comparison  with  the  structured  grid  solutions.  Details  of  the  grid  are  presented  in 
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Figure  24:  Pitching  airfoil. 

Table  3.  Based  on  a  temporal  refinement  study,  a  non-dimensional  time  step  [dt^  = 
tUoolc)  of  2.6  X  10“®  was  selected  for  the  computation. 

The  seven  cases  were  chosen  to  study  and  understand  the  effects  of  compressibility, 
pitch  rate  and  Reynolds  number  on  the  incipient  separation  process.  In  cases  1  through 
3  and  4  through  6,  the  pitch  rate  is  decreased  while  keeping  the  Reynolds  and  Mach 
numbers  constant.  Similarly  in  cases  1  and  4,  2  and  5,  3  and  6,  the  Mach  number  is 
increased  from  0.2  to  0.5  while  the  Reynolds  number  and  pitch  rate  are  fixed.  Cases  4 
and  7  are  computed  at  different  Reynolds  number  while  the  Mach  number  and  the  pitch 
rate  are  fixed.  These  computations  help  in  gaining  an  insight  into  several  important 
trends  related  to  the  change  in  Mach  number,  Reynolds  number  or  the  pitch  rate. 

8.2  Critical  Point  Theory 

The  understanding  of  unsteady  separation  for  a  pitching  airfoil  is  greatly  facilitated 
by  visualization  of  the  instantaneous  streamlines^.  There  are  several  possible  frames  of 

^Of  course,  the  instantaneous  streamlines  are  not  identical  to  the  particle  pathlines  or  streaklines 
since  the  flow  is  unsteady. 
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Table  3:  Details  of  the  grids 


Structured  Grid 


Case 

number 


Reynolds 

number 

Mach 

number 

Pitch 

rate 

10,000 

0.2 

0.2 

10,000 

0.2 

0.1 

10,000 

0.2 

0.05 

10,000 

0.5 

0.2 

10,000 

0.5 

0.1 

10,000 

0.5 

0.05 

100,000 

0.5 

0.2 

Ni 

Nr, 

Anfc 

xlO^ 

Asfc 

xlO^ 

Nbl 

637 

181 

1.0 

1.69 

64 

1011 

181 

1.0 

0.84 

64 

637 

325 

0.5 

1.69 

112 

1037 

181 

1.0 

1.025 

58 

2073 

361 

0.5 

0.511 

115 

1037 

181 

1.0 

1.025 

58 

2073 

361 

0.5 

0.511 

115 

1037 

181 

1.0 

1.025 

58 

2073 

361 

0.5 

0.511 

115 

1037 

181 

1.0 

1.025 

58 

2073 

361 

0.5 

0.511 

115 

2073 

159 

2.0 

0.509 

44 

4145 

317 

1.0 

0.252 

88 

2073 

337 

0.05 

0.509 

124 

4145 

673 

0.025 

0.252 

247 

Unstructured  Grid 

Case 

number 

Reynolds 

number 

Mach 

number 

Pitch 

rate 

Grid 

icell 

jnode 

Anjc 

xl0'‘ 

Asfc 

xlO^ 

Nbl 

1— 1 

10,000 

0.2 

0.2 

a 

30,238 

15,299 

10.0 

3.5 

14 

LEGEND  : 

= 

Nr,  = 
An  = 
As  = 
Nbl  = 
icell  = 
jnode  = 


Number  of  points  in  ^—direction 
Number  of  points  in  ?/— direction 
Average  normal  distance  of  points  next  to  airfoil 
Minimum  tangential  distance  of  points  on  airfoil 

Points  in  BL  measured  normal  to  airfoil  surface  at  mid-chord  for  a  =  0° 
Number  of  triangles 
Number  of  nodes 
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reference  including  1)  the  inertial  (laboratory)  frame  of  reference  wherein  the  freestream 
velocity  is  fixed  and  the  velocity  at  the  airfoil  surface  is  generally  nonzero,  and  2)  the 
rotating  frame  of  reference  (attached  to  the  airfoil)  wherein  the  freestream  velocity  is 
unsteady  and  the  velocity  at  the  airfoil  surface  is  zero.  The  rotating  frame  possesses 
the  seemingly  intuitive  advantage  for  physical  interpretation  of  “forward”  and  “reverse” 
flow  relative  to  the  airfoil.  Since  the  velocity  is  zero  at  the  airfoil  surface  in  the  rotating 
frame  of  reference,  the  fluid  immediately  adjacent  to  a  point  on  the  airfoil  surface  is 
either  instantaneously  moving  “forward”  (i.e.,  towards  the  trailing  edge)  or  “reverse” 
(i.e.,  towards  the  leading  edge). 

The  instantaneous  streamlines  in  the  rotating  frame  of  reference  were  selected  for 
analysis.  The  components  of  the  velocity  in  the  rotating  frame  of  reference  {x\y') 

are  related  to  the  components  of  the  velocity  (u,v)  in  the  inertial  frame  by  (Fig.  (24)) 

u'  =  p'  -f-  u  cos  0  +  n  sin  ^ 

v'  =  — fia:' —  Msin^  +  n  cos^  (99) 


Note  that  the  pitch  rate  is  negative  for  clockwise  (pitch  up)  rotation.  The  instanta¬ 
neous  streamlines  are 


^  =  u'(x',p',t) 

^  -  v'(x',y',t)  (100) 


where  r  is  the  parametric  length  along  the  instantaneous  streamlines  at  a  particular 
time. 


Equation  (100)  is  an  autonomous  system  of  ordinary  differential  equations  which  may 
possess  critical  points,  i.e.,  loci  where  u'  =  v'  —  0.  The  behavior  of  the  equations  in  the 
vicinity  of  these  critical  points  is  well  known  as  described,  for  example,  in  Kaplan  [17] 
and  Pontryagin  [31].  Recent  papers  by  Perry  and  Chong  [30]  and  Chong  et  al  [10]  have 
elucidated  further  details  regarding  fluid  motions  near  critical  points.  The  critical  points 


may  be  classified  on  the  basis  of  the  Jacobian  J  and  dilatation  A  where 

du'  dv'  du'  dv' 
dx'  dy'  dy'  dx' 

A  -  ^  ^ 

dx'  dy' 


(101) 

(102) 
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The  taxonomy  of  critical  points  for  two-dimensional  flow  is  described  in  Fig.  (25). 
For  J  <  0,  the  topology  is  a  saddle,  with  the  relative  orientation  of  the  asymptotes 
determined  by  the  value  of  A.  For  J  >  0  the  principal  topologies  are  nodes  and  foci, 
with  the  latter  obtained  for  J  >  A^/4.  The  patterns  are  stable  {i.e.,  the  instantaneous 
fluid  motion  is  towards  the  critical  point)  or  unstable  {i.e.,  the  instantaneous  fluid  motion 
is  away  from  the  critical  point)  based  on  the  sign  of  A.  Positive  A  leads  to  unstable 
patterns,  while  negative  values  of  A  lead  to  stable  patterns.  In  those  instances  where 
A  =  0,  the  topology  changes  from  a  saddle  (for  J  <  0),  to  a  simple  shear  (for  J  =  0) 
and  to  a  center  (for  J  >  0).  In  2-D  incompressible  flows,  the  dilatation  (A)  is  zero  and 
only  saddles  and  centers  can  occur.  Centers  appear  in  the  region  of  positive  Jacobian 
and  saddles  in  the  region  of  negative  Jacobian. 


y  =  J{(u,v)/(x,y)} 


Figure  25:  Classiflcation  of  critical  points. 

The  appearances  of  the  critical  points  is  closely  related  to  the  appearance  of  the 
recirculating  regions  and  the  development  of  the  flow  field.  In  the  following  section, 
critical  point  theory  is  applied  to  the  computation  of  unsteady  separation  on  a  pitching 
airfoil. 
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9  RESULTS 


The  flow  structures  observed  in  the  seven  cases  studied  are  described  in  the  present 
section.  A  linear  stability  analysis  was  also  conducted  as  a  part  of  the  research  program 
to  understand  the  reason  for  the  formation  of  the  recirculating  regions  in  the  flow  fleld. 
The  stability  analysis  has  been  presented  at  the  end  of  this  section. 

The  instantaneous  streamlines  obtained  from  the  computations  are  presented  in  the 
following  subsections.  The  results  presented  here  correspond  to  Grid  b  for  each  of  the 
cases  (except  for  the  unstructured  grid  results,  which  are  for  Grid  a).  The  instanta¬ 
neous  streamlines  are  plotted  based  on  a  reference  frame  attached  to  the  airfoil.  (The 
instantaneous  streamlines  have  not  been  plotted  to  scale  in  order  to  show  all  flow  de¬ 
tails.)  The  origin  of  the  coordinate  axis  is  the  quarter  chord  point.  The  airfoil-attached 
reference  frame  allows  an  unambiguous  definition  of  forward  and  reversed  flow  within  a 
thin  unseparated  boundary  layer,  since  the  velocity  of  the  fluid  at  the  airfoil  surface  is 
zero.  Forward  flow  is  defined  as  fluid  moving  toward  the  trailing  edge,  and  reversed  flow 
indicates  fluid  moving  toward  the  leading  edge. 

9.1  Case  1  :  Rcc  =  10“^,  Moo  =  0.2,  Q+  =  0.2  (Structured  Grid) 

The  instantaneous  streamlines  at  a  =  14.5°,  19.5°,  and  22.5°  are  shown  in  figures 
26  through  28.  At  a  =  14.5°,  the  flow  on  the  upper  surface  in  the  vicinity  of  the 
leading  edge  has  a  thin  reversed  flow  region  extending  to  the  7%  chord  position.  A 
clockwise-rotating  primary  recirculating  region  can  be  seen  near  the  leading  edge  on 
the  upper  surface  at  a  =  19.5°.  The  counterclockAvise- rotating  secondary  recirculating 
region  appears  beneath  the  primary  recirculating  region  and  a  clockwise-rotating  tertiary 
recirculating  region  between  the  leading  edge  and  the  primary  recirculating  region  at 
a  =  22.5°.  The  emergence  of  the  primary  recirculating  region  is  traced  to  the  appearance 
of  a  pair  of  critical  points  (which  are  defined  as  points  in  the  velocity  field  where  u=v=0 
in  a  reference  frame  attached  to  the  airfoil)  within  the  flowfield  at  a  =  14.99°  at  the 
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Figure  26:  Instantaneous  streamlines  at  a  =  14.5°  {Re^  =  10^,  =  0.2,  =  0.2; 

Case  1) 

18%  chord  location  and  a  distance  2.6  x  10“®c  above  the  airfoil  surface.  In  figure  29, 
the  close-up  of  the  instantaneous  streamlines  at  a  =  15°  show  the  structure  of  the 
primary  recirculating  region  just  after  its  birth.  The  instantaneous  streamlines  in  the 
region  between  forward  and  reversed  flow  are  in  the  shape  of  a  “C”.  As  a  increases,  one 
of  the  “C”  shaped  instantaneous  streamlines  collapses  at  a  single  point.  The  critical 
point  initially  appears  as  a  single  point  corresponding  to  a  pure  shear,  i.e.,  J  =  A  =  0 
(Fig.  25).  It  immediately  bifurcates  into  two  critical  points  (a  center  and  a  saddle) 
which  move  apart  and  the  center  gives  rise  to  the  primary  recirculating  region.  The 
primary  and  secondary  recirculating  regions  interact  with  each  other  to  eject  the  fluid 
close  to  the  airfoil  surface  in  a  direction  approximately  normal  to  the  wall,  thus  leading 
to  boundary  layer  separation  (Peridier  et  al  [28,  29]).  The  primary  recirculating  region 
eventually  detaches  from  the  airfoil  surface  and  becomes  the  dynamic  stall  vortex.  The 
developments  near  the  leading  edge  are  relatively  isolated  from  the  developments  on  the 
aft  portion  of  the  airfoil  through  a  region  of  attached  flow  (figure  30). 
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Figure  31:  Instantaneous  streamlines  at  a  =  14.5°  employing  unstructured  grid  {Rcc  = 
10^  Moc  =  0.2,  12+  =  0.2;  Case  1) 


9.2  Case  1  :  Rcc  =  10'^,  Moo  =  0.2,  17^  =  0.2  (Unstructured  Grid) 

Instantaneous  streamlines  obtained  using  the  unstructured  grid  code  are  shown  in 
Figs.  31  to  35  at  a  =  14.5°,  16.5°,  19.5°,  21.0°  and  22.5°.  Overall,  there  is  close  agree¬ 
ment  between  the  structured  and  unstructured  grid  computations.  The  primary  recir¬ 
culating  region  is  evident  in  both  cases  at  a  =  16.5°,  and  has  grown  at  a  =  19.5°.  The 
secondary  and  tertiary  recirculating  regions  at  a  =  22.5°  are  quantitatively  very  similar 
in  both  computations.  The  close  comparison  between  the  separate  computations  using 
two  fundamentally  different  algorithms  establishes  the  credibility  of  the  results. 
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a  =16.5” 


Figure  32:  Instantaneous  streamlines  at  «  =  16.5°  employing  unstructured  grid  [Rtc 
10^  Moo  ==  0.2,  =  0.2;  Case  1) 


a  =19.5° 


Figure  33:  Instantaneous  streamlines  at  a  =  19.5°  employing  unstructured  grid  {Rec 
10^  Moo  =  0.2,  =  0.2;  Case  1) 
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a  =  21.0° 


Figure  34:  Instantaneous  streamlines  at  a  =  21.0°  employing  unstructured  grid  {Rec 
10^  Moo  =  0.2,  =  0.2;  Case  1) 


a  =  22.5° 


Figure  35:  Instantaneous  streamlines  at  a  =  22.5°  employing  unstructured  grid  {Rcc 
10^  Moo  =  0.2,  =  0.2;  Case  1) 
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9.3  Case  2  :  Rcc  =  lOS  =  0.2,  Q+  =  0.1 

The  instantaneous  streamlines  are  shown  in  figures  36  through  39  for  a  =  13.5°,  15°, 
17°  and  18°,  respectively.  At  a  =  13.5°,  there  is  reversed  flow  on  the  upper  surface 
extending  almost  to  the  leading  edge.  However,  there  are  no  recirculating  regions  in 
the  area  shown.  (Flow  near  the  trailing  edge  shows  complex  separation  behavior.)  The 
clockwise  rotating  primary  recirculating  region  forms  at  ct  =  14.29°  at  27%  chord  and 
a  distance  1.25  x  lO^^c  above  the  wall.  At  a  =  15°,  the  primary  recirculating  region 
is  shown.  At  a  =  17°,  a  counterclockwise  rotating  secondary  recirculating  region  forms 
below  the  primary  recirculating  region  at  the  16%  chord  position.  The  primary  recircu¬ 
lating  region  has  grown  normal  to  the  airfoil  surface  and  also  moved  toward  the  leading 
edge.  The  instantaneous  streamlines  are  shown  over  the  entire  airfoil  at  a  =  18°  in 
figure  39,  which  shows  that  the  developments  in  the  leading  edge  region  is  relatively 
isolated  from  developments  in  the  trailing  edge  region.  A  clockwise  rotating  tertiary 
recirculating  region  develops  on  top  of  the  secondary  recirculating  region  in  between  the 
leading  edge  and  the  primary  recirculating  region.  The  flow  structures  are  similar  to 
Case  1  (fi5+  =  0.2).  However,  the  primary  recirculating  region  appears  at  a  smaller  angle 
as  compared  to  Case  1.  Also,  the  structures  form  farther  away  from  the  leading  edge 
and  toward  the  trailing  edge  on  the  upper  surface  in  Case  2. 
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Figure  36:  Instantaneous  streamlines  at  a  =  13.5°  {Rbc  =  in'*,  Moo  =  0.2,  =  0.1; 

Case  2) 


Figure  37:  Instantaneous  streamlines  at  a  =  15°  {Rcc  =  10^,  Moo  =  0.2,  12+  =  0.1;  Case 
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Figure  38:  Instantaneous  streamlines  at  a  =  17°  {Rtc  =  10^,  Moo  =  0.2,  11+  =  0.1;  Case 
2) 
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Figure  39:  Instantaneous  streamlines  over  the  entire  airfoil  at  a  =  18°  [Rtc  =  10“*, 
Moo  =  0.2,  n+  =  0.1;  Case  2) 


9.4  Case  3  :  Rec  =  lOS  Moo  =  0.2,  Q+  =  0.05 

The  instantaneous  streamlines  over  the  entire  airfoil  are  shown  in  figures  40  through 
43  for  a  =  9°,  10.5°,  12°,  and  13.5°,  respectively.  At  a  =  9°,  the  reversed  flow  extends 
over  most  of  the  airfoil  upper  surface  and  a  trailing  edge  recirculating  region  can  be  seen. 
The  reversed  flow  extends  over  the  entire  airfoil  upper  surface  at  a  =  10.5°.  The  trailing 
edge  recirculating  region  stretches  along  the  airfoil  upper  surface  and  breaks  down  into 
multiple  recirculating  regions  at  a  =  12°  extending  over  the  entire  airfoil  top  surface. 
Unlike  the  Cases  1  and  2  (which  are  at  higher  pitch  rate),  the  flow  developments  near 
the  trailing  edge  are  not  isolated  from  the  developments  near  the  leading  edge.  Also, 
the  recirculating  regions  appear  far  from  the  wall  (average  distance  of  4.0  x  10~^c).  The 
instantaneous  streamlines  at  a  =  13.5°  demonstrate  the  increase  in  size  of  the  multiple 
recirculating  regions  in  the  transverse  direction  with  increase  in  the  angle  of  attack, 
giving  rise  to  breakdown  of  the  entire  upper  surface  boundary  layer.  The  development 
of  the  flow  structure  which  gives  rise  to  the  dynamic  stall  vortex  is  not  isolated  in  the 
leading  edge  only  and  appears  over  the  entire  airfoil  upper  surface,  which  is  different 
from  Cases  1  and  2. 
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Figure  44:  Instantaneous  streamlines  at  a  =  18°  {Rec  =  10^,  Moo  =  0.5,  =  0.2;  Case 

4) 

9.5  Case  4  :  Rcc  =  104,  Moo  =  0.5,  —  0-2 

The  instantaneous  streamlines  over  the  front  50%  of  the  airfoil  are  shown  in  figures 
44  through  47  for  a  =  18°,  19.5°,  22.5°,  and  25.5°,  respectively.  The  reversed  flow 
region  observed  over  the  airfoil  upper  surface  at  a  =  18°  is  much  thicker  as  compared 
to  the  reversed  flow  region  at  Moo  =  0.2  (Cases  1  and  2).  At  a  =  19.5°,  the  primary 
recirculating  region  has  formed.  In  particular,  the  primary  recirculating  region  appears 
first  at  a  =  18.8°  at  30%  chord  position  and  a  distance  1.7  x  10“^c  from  the  airfoil  surface. 
The  secondary  recirculating  region  can  be  seen  at  a  =  22.5°.  The  secondary  recirculating 
region  develops  below  the  primary  recirculating  region  which  expands  normal  to  airfoil 
surface.  At  a  =  25.5°,  the  tertiary  recirculating  region  appears  above  the  secondary 
recirculating  region  and  the  primary  recirculating  region  is  pushed  farther  from  the 
airfoil  surface.  In  this  case,  the  vortical  structures  have  larger  length  scale  as  compared 
to  Case  1  or  Case  2.  Again,  the  development  of  the  flow  field  is  very  similar  to  Cases  1 
and  2,  but  with  a  delay. 
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Figure  45;  Instantaneous  streamlines  at  a  =  19.5°  [Rcc  =  10^,  il 
Case  4) 


Figure  47:  Instantaneous  streamlines  at  a  =  25.5°  [Rcc  =  10^,  Moo  —  0.5,  =  0.2; 

Case  4) 

9.6  Case  5  :  =  lOS  Moo  =  0.5,  Q+  =  0.1 

The  instantaneous  streamlines  at  a  =  15°,  16.5°,  18°,  and  19.5°  over  the  leading  70% 
of  the  airfoil  are  shown  in  figures  48  through  51.  At  a  =  15°,  a  thick  reversed  flow 
region  is  observed  over  the  entire  upper  surface.  The  primary  recirculating  originates  at 
a  =  15.65°  at  44%  chord  position  and  a  distance  3.8  x  10“^c  above  the  airfoil  surface. 
The  primary  recirculating  region  is  observed  to  have  developed  at  a  =  16.5°,  and  the 
secondary  recirculating  region  forms  below  the  primary  recirculating  region  (a  =  18°). 
At  a  =  19.5°,  the  tertiary  recirculating  region  can  be  seen.  The  secondary  recirculating 
region  separates  the  primary  recirculating  region  from  the  tertiary  recirculating  region. 
The  separation  process  is  again  observed  to  be  similar  to  Cases  1,  2,  and  4.  The  principal 
differences  are  that  the  structures  form  farther  away  from  the  leading  edge  compared 
to  the  other  three  cases,  and  the  length  scale  associated  with  the  flow  structures  is 
substantially  greater. 
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Figure  50:  Instantaneous  streamlines  at  a  =  18°  (Rec  =  10^,  M^o  =  0.5,  =  0.1;  Case 
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Figure  51:  Instantaneous  streamlines  at  a  =  19.5°  {Rcc  =  10^,  Moo  =  0.5,  =  0.1: 
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Figure  52;  Instantaneous  streamlines  at  a  =  10.5°  {Re^  —  10^,  Moo  =  0-5,  =  0.05; 

Case  6) 

9.7  Case  6  :  Rec  =  lOS  Moo  =  0.5,  Q+  =  0.05 

The  instantaneous  streamlines  over  the  entire  airfoil  upper  surface  are  shown  in  fig¬ 
ures  52  through  55  for  a  =  10.5°,  12°,  13.5°,  and  15°,  respectively.  At  a  =  10.5°, 
a  recirculating  region  is  present  near  the  trailing  edge  and  the  reversed  flow  extends 
over  the  entire  airfoil  upper  surface.  The  trailing  edge  recirculating  region  divides  into 
two  separate  recirculating  regions  at  a  =  12°.  The  recirculating  region  attached  at  the 
trailing  edge  is  shed  in  the  wake  at  a  =  13.5°,  while  the  other  one  expands  and  is  still 
attached  to  the  airfoil  surface.  With  increase  in  a  two  more  recirculating  regions  appear 
and  are  observed  at  a  =  15°.  The  development  of  the  flow  structures  in  this  case  resem¬ 
bles  Cases  1,  2,  4,  and  5.  However,  the  flow  structures  are  not  constrained  to  only  a  part 
of  the  airfoil  upper  surface.  The  primary  recirculating  region  appears  near  the  trailing 
edge  and  extends  over  the  entire  airfoil  surface  before  breaking  down  into  secondary  and 
tertiary  recirculating  regions. 
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Figure  53:  I 
Case  6) 


Figure  54: 
Case  6) 


Figure  55:  Instantaneous  streamlines  at  o;  =  15°  [Re^  =  10^,  M^o  —  0.5,  12+  =  0.05; 
Case  6) 

9.8  Case  7  :  Rcc  =  10^,  Mqo  =  0.5,  QJ  =  0.2 

The  instantaneous  streamlines  are  presented  for  a  =  15°,  18°,  19°,  and  19.5°  in 
figures  56  through  59.  The  instantaneous  streamlines  at  ct  =  15°  (figure  56)  show  the 
presence  of  the  primary  recirculating  region  near  the  leading  edge  in  the  close-up  figure 
of  the  boxed  region.  In  particular,  the  primary  recirculating  region  forms  at  a  =  14.9° 
as  compared  to  a  =  18.8°  for  Case  4.  The  primary  recircidating  region  expands  in  a 
direction  normal  to  airfoil  surface  with  increase  in  a.  This  can  be  observed  at  a  =  18°. 
At  a  =  19°,  the  secondary  and  tertiary  recirculating  regions  are  also  evident.  Up  to 
a  =  19°,  the  development  of  the  flow  field  is  very  similar  to  Cases  1,  2,  4,  5  and  6. 
However,  at  a  =  19.5°,  midtiple  recirculating  regions  are  seen  attached  to  the  airfoil 
surface.  Similar  residts  have  been  reported  in  an  analytical  study  by  Smith  [38]  and  in  a 
numerical  study  by  Ghia  et  al  [15]  at  Rcc  =  4.5  x  10“*,  where  simultaneous  appearances 
of  multiple  recirculating  regions  near  the  leading  edge  at  high  Reynolds  numbers  were 
found.  The  pressure  coefficient  contours  at  a  =  18°,  and  19.5°  in  figures  60  and  61, 
respectively,  show  the  appearance  of  a  shock  in  the  flowfield  at  a  =  19.5°.  Examination 
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Figure  56:  Instantaneous  streamlines  at  a  =  15°  {Rcc  =  10^,  Moo  —  0.5,  11+  =  0.2;  Case 

7) 

of  the  Mach  contours  and  the  divergence  of  velocity  in  the  flowfield  (not  shown  here) 
also  indicates  a  shock.  The  shock  appears  between  a  =  19°  and  a  =  19.5°.  A  separate 
inviscid  computation  showed  similar  appearance  of  a  shock  near  the  leading  edge,  which 
proves  that  the  appearing  shock  is  an  inviscid  phenomenon.  For  the  inviscid  case,  the 
shock  appears  at  a  smaller  angle  of  attack  and  nearer  to  the  leading  edge  than  in  the 
viscous  case.  The  length  scale  associated  with  the  recirculating  regions  in  this  case  is 
observed  to  be  very  small  as  compared  to  those  in  the  other  cases. 

9.9  Effect  of  compressibility 

The  results  indicate  important  trends  related  to  the  increase  in  Mach  number.  In  the 
Cases  1  &  4  and  2  &  5,  the  pitch  rate  and  Reynolds  number  are  fixed  while  the  Mach 
number  is  increased  from  0.2  to  0.5.  Increasing  the  Mach  number  (at  fixed  Reynolds 
number  and  pitch  rate)  causes  the  primary  recirculating  region  to  form  farther  from 
the  airfoil  surface  and  delays  its  formation.  Sankar  and  Tassa  [35]  also  saw  a  delay 
in  the  formation  of  the  “leading  edge  vortex”  (primary  recirculating  region)  when  the 
Mach  number  was  increased  from  0.2  to  0.4  at  Rcc  =  5.0  x  10^.  Figure  62  displays 
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Figure  61:  Pressure  coefficient  contours  at  a=19.5°  {Rbc  =  10®,  M^o  —  0.5,  12+  =  0.2; 
Case  7) 


the  angle  at  which  the  primary  recirculating  region  first  appears  as  a  function  of  the 
pitch  rate  for  two  different  Mach  numbers  at  Rbc  =  10^.  At  a  constant  pitch  rate,  the 
formation  of  the  primary  recirculating  region  is  delayed  to  higher  angles  when  the  Mach 
number  is  increased  from  0.2  to  0.5.  The  surface  pressure  coefficient  was  compared 
to  study  the  effect  of  Mach  number.  Figure  63  shows  that  for  a  fixed  pitch  rate  and 
Reynolds  number,  increasing  the  Mach  number  causes  a  decrease  in  the  leading  edge 
suction  pressure  coefficient  which  in  turn  results  in  a  lower  adverse  pressure  gradient  on 
the  airfoil  upper  surface.  The  decreased  adverse  pressure  gradient  on  the  top  surface 
at  higher  Mach  numbers  retards  the  movement  of  the  reversed  flow  region  toward  the 
leading  edge  but  increases  the  thickness  of  the  reversed  flow  region  and  eventually  delays 
the  formation  of  the  primary  recirculating  region.  The  length  scale  associated  with  the 
flow  structures  increase  as  well.  The  decrease  in  the  magnitude  of  the  leading  edge 
suction  pressure  coefficient  is  primarily  due  to  the  compressibility  effects  on  the  boundary 
layer,  which  increase  the  effective  radius  of  the  leading  edge  of  the  airfoil  and  in  turn 
lead  to  a  decrease  in  the  leading  edge  suction  pressure  coefficient.  This  is  opposite  to 
the  behavior  expected  from  a  quasi-stationary  application  of  the  Prandtl-Glauert  rule. 
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Figure  62:  Effect  of  Mach  number  and  pitch  rate  on  the  appearance  of  the  primary 
recirculating  region  at  Rtc  —  10“^ 

where  Cpj^fCp-  =  l/\/l  —  pressure  coefficient  for  a  compressible  case  of 

Mach  number  M  and  Cp^  is  the  pressure  coefficient  for  an  incompressible  case). 

9.10  Effect  of  pitch-rate 

The  flowfield  is  strongly  dependent  on  the  pitch  rate.  Comparison  of  Cases  1,  2  and 
3,  or  4,  5  and  6  exhibit  the  effects  of  the  pitch  rate  at  a  fixed  Mach  number  and  Reynolds 
number.  The  pitch  rate  is  decreased  from  0.2  through  0.05  in  cases  1  through  3  and  4 
through  6.  Increase  in  the  pitch  rate  of  the  airfoil  delays  the  formation  of  the  primary 
recirculating  region  to  higher  angles  of  attack.  The  primary  recirculating  region  forms 
closer  to  the  leading  edge  on  the  upper  surface.  It  is  also  observed  that  decrease  in 
the  pitch  rate  causes  the  transition  from  a  predominantly  leading-edge  separation  to  a 
complex  separation  over  the  entire  airfoil.  Figure  62  shows  the  delay  in  the  formation 
of  the  recirculating  region  with  increase  in  pitch  rate  at  constant  Mach  numbers  and 
Rcc  —  10'*.  It  should  be  noted  that  the  delay  in  the  formation  of  the  recirculating  region 
to  higher  angles  does  not  mean  that  there  is  an  increase  in  the  actual  time  from  start 
of  pitch- up  to  the  appearance  of  the  recirculating  region  (the  pitch  rate  is  different). 
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Figure  63:  Comparison  of  the  surface  pressure  coefficient  at  a  =  16.5°  for  two  different 
Mach  numbers  {Rec  =  10^,  12+  =  0.2) 

At  higher  pitch  rates,  there  is  a  decrease  in  the  adverse  pressure  gradient  on  the  upper 
surface  of  the  airfoil  (figure  64)  that  can  also  be  explained  from  the  potential  flow  theory, 
which  slows  down  the  development  of  the  reversed  flow  region  and  defers  the  appearance 
of  the  primary  recirculating  region  to  a  higher  angle. 

9.11  Effect  of  Reynolds  number 

The  flowfleld  is  sensitive  to  the  change  in  Reynolds  number.  Cases  4  and  7  are 
compared  to  study  the  effect  of  increase  in  Reynolds  number  at  a  fixed  Mach  number 
and  pitch  rate.  Increase  in  the  Reynolds  number  hastens  the  appearance  of  the  primary 
recirculating  region  and  decreases  the  length  scale  of  the  flow  structures.  The  primary 
recirculating  region  forms  closer  to  the  leading  edge  on  the  upper  surface  at  higher 
Reynolds  number.  At  Rcc  —  10®  and  Moo  =  0.5,  a  shock  appears  on  the  top  surface.  The 
appearance  of  the  shock  was  found  to  be  an  inviscid  phenomenon.  Multiple  recirculating 
regions  develop  simultaneously  near  the  leading  edge  at  higher  Reynolds  number  and 
cause  a  complex  separation  of  the  boundary  layer.  Up  to  the  formation  of  the  shock,  the 
evolution  of  the  flowfleld  is  identical  in  both  cases  except  for  the  lag.  The  formation  of  the 
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Figure  64:  Comparison  of  the  surface  pressure  coefficient  at  a  =  13.5°  for  two  different 
pitch  rates  (iJco  =  10“*,  M^o  =  0.2) 

recirculating  regions  near  the  leading  edge  is  not  induced  by  the  shock.  Chandrasekhara 
et  al  [9]  observed  the  presence  of  similar  shocks  near  the  leading  edge  at  Moo  =  0.45  and 
Rcc  =  3.6  X  10®,  in  the  experiments. 

9.12  Linear  Stability  Analysis 

A  linear  incompressible  stability  analysis  was  conducted  to  determine  whether  the 
formation  of  the  primary  recirculating  region  is  related  to  an  instability  of  the  flow.  The 
basic  flow  required  for  the  stability  analysis  was  obtained  from  the  numerical  simula¬ 
tions  for  Cases  1  and  2,  which  have  been  described  in  the  preceding  subsections.  The 
Mach  number  in  the  flowfleld  does  not  exceed  0.6  in  both  the  cases,  and  therefore  an 
incompressible  analysis  was  deemed  sufficient. 

In  the  stability  analysis  of  a  basic  flow  field,  a  small  disturbance  is  introduced  to  the 
basic  flow  and  the  development  of  the  flow  field  in  time  is  sought  to  determine  whether 
the  flow  field  is  stable  or  unstable.  The  velocities  and  the  pressure  are  written  as 

u  =  U-\-u'\  v  —  V^v'\  p  =  P p'  (103) 

where  U,  V,  and  P  denote  the  basic  flow  and  u',  and  p'  denote  the  disturbances  in  the 
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velocities  and  the  pressure.  These  expressions  are  substituted  into  the  incompressible 
Navier-Stokes  equations  which  are  then  linearized.  The  disturbances  are  described  as  a 
complex  perturbation  streamfunction  of  the  form 

y-,  i)  =  Hy)  exp{i{fca:  -  ujt)}  (104) 


where  k  and  u;  are  the  complex  wavenumber  and  frequency  respectively,  (p{y)  is  the  com¬ 
plex  amplitude  function.  The  perturbed  velocities  are  defined  in  terms  of  the  perturbed 
streamfunction  (■0')  and  the  pressure  is  eliminated  from  the  linearized  Navier-Stokes 
equation,  which  result  in  the  following  equation 
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-  -{D^  -  k'^)ip  =  0  (105) 


where  the  operator  ‘‘D'  denotes  differentiation  with  respect  to  y.  Details  of  the  derivation 
of  the  above  equation  is  given  in  Appendix  E. 

Equation  (105)  is  a  fourth  order  complex  differential  equation  which  together  with 
the  appropriate  boundary  condition  defines  the  stability  eigenvalue  problem.  In  the 
present  study,  the  reference  frame  is  assumed  to  be  attached  to  the  airfoil  surface  and 
the  X-  and  y-directions  are  aligned  in  the  tangential  and  normal  directions  to  the  airfoil 
surface,  respectively.  The  airfoil  is  pitching  up  at  a  constant  pitch  rate  and  therefore 
the  boundary  conditions  can  be  written  as 


u'  =  v'  =  Q  or  ’f{y)  =  D<p{y)  =  0  at  y  —  0 


(106) 


<Pu'  dv'  /  \  T-.'?  /  ^ 

=  ^  =  D^^{y)  =  0  at  y  =  oo  (107) 

The  stability  eigenvalue  problem  was  solved  by  utilizing  the  Chevyshev  polynomials 
to  approximate  the  eigenfunctions  [59].  To  use  the  Chevyshev  polynomial  as  a  basis, 
the  domain  [0,oo]  is  mapped  into  the  interval  [-1,1]  with  the  new  variable  z  defined  as 
2:  =  1  —  2exp{—y/yo),  where  yo  is  the  edge  of  the  boundary  layer.  The  eigenvalues  of  the 
system  of  equation  are  determined  using  the  IMSL  subroutines  and  the  most  unstable 
mode  selected  out  of  the  set  of  eigenvalues,  based  on  the  imaginary  part  of  the  wave 
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frequency  (a;)  with  the  highest  value.  The  flow  is  unstable  when  the  imaginary  part  of 
the  wave  frequency  is  positive  and  its  magnitude  is  a  measure  of  the  temporal  growth 
rate  of  the  disturbance. 

The  developed  solver  was  validated  by  performing  two  computations  to  ascertain 
the  accuracy.  Stability  computations  were  performed  to  determine  the  neutral  curve 
(lm[c<;]=0)  for  a  Blasius  profile  and  a  Poiseuille  flow.  Excellent  agreement  was  found 
with  previous  computations  for  both  the  cases. 

The  instability  of  the  flow  can  be  of  two  types  -  absolute  instability  and  convective 
instability.  A  flow  is  said  to  be  absolutely  unstable  if  localized  disturbances  spread 
upstream  and  downstream  and  influence  the  entire  flow  field.  On  the  other  hand,  if 
the  disturbances  are  swept  away  from  the  source  but  nevertheless  increase  in  magnitude 
with  time,  the  flow  is  said  to  be  convectively  unstable.  Figure  65  shows  schematically 
the  two  different  types  of  instabilities. 

To  study  the  absolute  instability,  complex  wavenumbers  are  selected  and  the  cor¬ 
responding  most  unstable  wave  frequencies  are  computed  to  determine  the  dispersion 
relation  of  complex  k  and  complex  a>.  The  point  in  the  k  —  u  plane  (wq,  K)  where 
dujfdk  =  0  gives  the  branch  point  singularity  of  the  dispersion  relation.  The  disturbance 
will  grow  or  decay  at  its  location  of  generation  depending  on  the  sign  of  Imfwo]-  The 
flow  is  absolutely  unstable  if  the  branch  point  singularity  of  its  dispersion  relation  lies 
in  the  upper  half  of  the  complex  oj  plane  (Im[c4;o]>0).  To  study  convective  instability  of 
the  flow,  only  real  values  of  the  wavenumbers  are  selected  and  the  most  unstable  wave 
frequency  computed.  The  disturbance  wiU  grow  in  time  at  a  location  away  from  the 
source  (convectively  unstable)  if  lm[c(;]>0. 

The  stability  of  the  flowfield  was  computed  and  analyzed.  The  velocity  profiles  at  dif¬ 
ferent  chordwise  locations  and  angles  of  attack,  obtained  from  the  numerical  simulation, 
were  studied.  In  Case  1,  the  primary  recirculating  region  appears  at  approximately  15°. 
The  stability  analysis  for  this  case  showed  that  the  flowfield  is  not  absolutely  unstable  at 
any  location  in  the  front  50%  of  the  chord  up  to  an  angle  of  attack  of  18°,  which  is  well 
past  the  angle  at  which  the  primary  recirculating  region  appears.  This  demonstrates 
that  the  primary  recirculating  region  does  not  form  due  to  an  absolute  instability  of  the 
flow. 
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Figure  65:  Sketch  of  a  typical  impulse  response  :  (a)  stable;  (b)  absolutely  unstable;  (c) 
convectively  unstable 
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Figure  66:  Majdmum  temporal  growth  rate  of  the  disturbance  in  the  flowfield  at  various 
chordwise  positions  on  the  airfoil  and  at  different  angles  of  attack  (Case  1) 

Convective  stability  analysis  was  conducted  for  Case  1.  Figure  66  shows  the  maxi¬ 
mum  temporal  growth  rate  for  the  velocity  profiles  at  different  chordwise  positions  of  the 
airfoil  for  different  angles  of  attack.  Figure  67  shows  the  wavenumbers  corresponding  to 
the  maximum  temporal  growth  rate  at  different  chord  positions  of  the  airfoil.  It  is  seen 
that  the  flow  is  convectively  unstable  over  the  entire  airfoil  upper  surface  and  the  insta¬ 
bility  increases  in  magnitude  with  the  increase  in  the  angle  of  attack.  The  magnitude  of 
convective  instability  is  also  higher  in  the  region  near  the  leading  edge  and  increases  at 
a  faster  rate  with  increase  in  the  angle  of  attack  as  compared  to  the  aft  portion  of  the 
airfoil.  From  the  numerical  simulation  we  had  seen  that  the  primary  recirculating  region 
appears  around  18%  chord  position  for  Case  1.  Therefore  it  can  be  said  that  the  appear¬ 
ance  of  the  primary  recirculating  region  may  be  due  to  the  convective  instability  of  the 
flowfield  near  the  leading  edge.  The  wavenumber  of  the  disturbance  (normalized  by  the 
chord  length)  corresponding  to  the  maximum  temporal  growth  rate  at  the  18%  chord 
position  is  approximately  50  (figure  67),  which  corresponds  to  a  wavelength  of  around 
12.5%  of  the  chord  and  corresponds  roughly  to  the  size  of  the  primary  recirculating 
region. 

Figures  68  and  69  show  the  maximum  temporal  growth  rate  and  the  corresponding 
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Figure  67:  Wavenumber  corresponding  to  maximum  temporal  growth  rate  of  the  dis¬ 
turbance  at  various  chordwise  positions  on  the  airfoil  and  at  different  angles  of  attack 
(Case  1) 

wave  number  at  different  chordwise  positions  and  angles  of  attack  for  Case  2.  Case  2 
represents  a  decrease  in  the  pitch  rate  as  compared  to  Case  1  while  the  Mach  number 
and  Reynolds  number  are  the  same.  It  is  observed  that  the  magnitude  of  the  instability 
(temporal  amplification  rate)  is  much  higher  for  Case  2  as  compared  to  Case  1  at  any 
particular  angle.  This  suggests  that  the  primary  recirculating  region  will  appear  at  a 
smaller  angle  of  attack  for  this  case  with  respect  to  Case  1.  A  similar  trend  was  observed 
in  the  numerical  simulations  which  is  shown  in  figure  62.  The  wavenumber  (normalized 
by  the  chord  length)  corresponding  to  the  maximum  growth  rate  at  the  location  where 
the  primary  recirculating  region  appears  is  approximately  40,  which  in  turn  relates  to  a 
disturbance  wavelength  of  around  16%  of  the  chord.  Sections  9.1  and  9.3  describe  that 
the  decrease  in  the  pitch  rate  at  a  constant  Mach  number  and  Reynolds  number  leads 
to  an  increase  in  the  size  of  the  primary  recirculating  region,  which  is  also  indicated  in 
the  linear  stability  analysis. 

The  linear  stability  analysis  shows  that  the  appearance  of  the  primary  recirculating 
region  is  related  to  the  convective  instability  of  the  flowfield  near  the  leading  edge  of 
the  airfoil.  The  wavelength  of  the  disturbance  is  comparable  to  the  size  of  the  primary 
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Figure  68:  Mciximum  temporal  growth  rate  of  the  disturbance  in  the  flowfield  at  various 
chordwise  positions  on  the  airfoil  and  at  different  angles  of  attack  (Case  2) 


Figure  69:  Wavenumber  corresponding  to  maximum  temporal  growth  rate  of  the  dis¬ 
turbance  at  various  chordwise  positions  on  the  airfoil  and  at  different  angles  of  attack 
(case  2) 
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recirculating  region  as  seen  from  the  full  Navier-Stokes  simulations.  The  trend  with 
regards  to  the  effect  of  change  in  pitch  rate  on  the  size  of  the  primary  recirculating 
region  and  the  angle  at  which  it  first  appears  matches  with  the  observations  from  the 
direct  numerical  simulation. 

9.13  Flow  Control 

Control  of  dynamic  stall  is  very  important  to  effectively  utilize  it  in  unsteady  aerody¬ 
namic  applications.  There  have  been  very  few  studies  in  the  past  to  investigate  different 
flow  control  techniques  to  delay  dynamic  stall.  The  active  control  using  modulated  suc¬ 
tion/injection  to  manage  the  dynamic  stall  vortex  has  been  studied  the  most  and  has 
been  quite  successful.  But,  there  are  practical  problems  in  incorporating  this  control 
technique  in  the  small  blades  of  the  helicopter  rotor. 

The  present  research  project  has  been  successful  in  tracing  the  formation  of  the 
dynamic  stall  vortex  to  a  critical  point  in  the  flowfield.  This  information  may  be  very 
useful  in  investigating  different  flow  control  methods  because  it  gives  an  idea  with  respect 
to  the  location  in  the  flowfield  where  flow  control  has  to  be  applied  to  be  most  effective. 

A  study  was  conducted  to  delay  the  formation  of  the  primary  recirculating  region 
by  preferential  heating  and  cooling  of  the  entire  airfoil  surface.  It  has  been  observed  in 
the  results  presented  that  the  formation  of  the  primary  recirculating  region  depends  to 
a  great  extent  on  the  magnitude  of  the  adverse  pressure  gradient  on  the  airfoil  upper 
surface.  At  higher  Mach  numbers,  the  leading  edge  suction  pressure  coefficient  decreases 
as  compared  to  that  at  lower  Mach  numbers.  This  leads  to  a  decrease  in  the  adverse 
pressure  gradient  on  the  airfoil  upper  surface  and  an  eventual  delay  in  the  formation  of 
the  primary  recirculating  region.  The  same  idea  of  decreasing  the  leading  edge  suction 
pressure  coefficient  can  be  utilized  to  delay  the  formation  of  the  primary  recirculating 
region.  A  preferential  heating  of  the  region  near  the  leading  edge  and  cooling  of  the 
trailing  edge  region  can  effectively  lead  to  a  decrease  in  the  leading  edge  pressure  coeffi¬ 
cient  and  increase  in  pressure  coefficient  near  the  trailing  edge.  This  on  the  other  hand 
may  lead  to  a  decrease  in  the  adverse  pressure  coefficient  on  the  airfoil  upper  surface. 

Simultaneous  heating  of  the  leading  edge  region  and  cooling  of  the  trailing  edge 
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Figure  70:  Profile  of  the  temperature  applied  on  airfoil  surface 

region  was  applied  to  Case  1  {Rec  —  10^,  Moo  =  0.2,  =  0.2).  Both  upper  and 

lower  surfaces  of  the  airfoil  were  either  heated  or  cooled  for  simplicity.  The  temperature 
profile  that  was  applied  on  the  airfoil  surface  is  shown  as  a  function  of  the  chordwise 
position  in  figure  70.  The  profile  is  a  modified  sinusoidal  function.  The  temperature 
was  applied  to  the  surface  at  the  start  of  the  pitch  up  motion  of  the  airfoil  from  zero 
degrees  as  an  exponential  function  of  time  to  reach  the  final  value  shown  in  the  figure 
at  approximately  1°. 

The  instantaneous  streamlines  over  the  airfoil  for  Case  1  with  modulated  heating 
and  cooling  are  shown  in  figures  71  to  75  at  a  =  13.5°,  15.0°,  16.5°,  18.0°,  and  21.0°.  At 
a  =  13.5°,  there  is  already  reversed  flow  over  the  entire  airfoil  upper  surface.  There  are 
also  two  recirculating  regions  present  near  the  trailing  edge  but  none  near  the  leading 
edge.  By  a  =  15.0°,  the  reversed  flow  has  considerably  increased  in  thickness  and  is  of 
the  order  of  the  thickness  of  the  airfoil.  But,  there  is  no  presence  of  primary  recirculating 
region  near  the  leading  edge.  For  Case  1  without  heating,  it  was  seen  that  the  primary 
recirculating  region  appeared  at  a  =  15.0°  at  18%  chord  position  and  the  reversed  flow 
was  much  thinner.  Figure  73  shows  the  presence  of  the  primary  recirculating  region  near 
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Figure  71:  Instantaneous  streamlines  at  a  =  13.5°  for  Case  1  with  modulated  heating 
and  cooling 

the  leading  edge  at  a  =  16.5°.  In  fact,  the  primary  recirculating  region  appears  due  to 
a  breakup  of  the  trailing  edge  vortex  and  at  approximately  60%  chord  position.  The 
primary  recirculating  region  grows  in  size  with  increase  in  the  angle  of  attack  and  moves 
upstream.  This  is  shown  in  figure  74  at  a  =  18.0°.  The  secondary  recirculating  region 
can  be  seen  at  a  =  21.90°  in  figure  75.  The  development  of  the  flowfield  is  very  similar 
to  the  case  where  there  is  no  modulated  heating  or  cooling.  The  difference  in  the  two 
cases  is  the  position  of  the  primary  recirculating  region  and  the  angle  of  attack  of  its 
first  appearance.  Modulated  heating  delays  the  formation  of  the  primary  recirculating 
region  by  about  1°.  The  recirculating  region  also  appears  at  a  farther  distance  from  the 
leading  edge  and  the  airfoil  upper  surface  as  compared  to  the  case  without  any  heating. 
The  flow  structures  are  much  bigger  in  size  and  cover  the  entire  airfoil  upper  surface. 

The  surface  pressure  coefficients  are  shown  in  figures  76  to  79.  The  figures  show 
a  comparison  of  the  surface  pressure  coefficients  for  Case  1  without  heating,  Case  1 
with  modulated  heating  and  cooling,  and  Case  4  (which  denotes  the  effect  of  increase 
in  compressibility).  From  the  figures  it  can  be  seen  that  modulated  heating  of  the 
leading  edge  region  and  cooling  of  the  trailing  edge  region  of  the  airfoil  surface  leads  to 
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Figure  72:  Instantaneous  streamlines  at  a  =  15.0°  for  Case  1  with  modulated  heating 
and  cooling 


Figure  73:  Instantaneous  streamlines  at  a  =  16.5°  for  Case  1  with  modulated  heating 
and  cooling 


Figure  74:  Instantaneous  streamlines  at  cc  =  18.0°  for  Case  1  with  modulated  heating 
and  cooling 


Figure  75:  Instantaneous  streamlines  at  ct  =  21.0°  for  Case  1  with  modidated  heating 
and  cooling 


Figure  76:  Surface  pressure  coefficient  at  a  =  13.5° 

a  substantial  drop  in  the  leading  edge  suction  coefficient  which  is  comparable  to  that 
at  Moo  =  0.5.  The  adverse  pressure  gradient  on  the  airfoil  surface  is  also  less  when  the 
airfoil  surface  is  heated.  But,  as  is  seen  from  the  instantaneous  streamline  plots,  this 
does  not  prevent  the  formation  of  the  reversed  flow  region  on  the  airfoil  upper  surface. 

To  summarize,  it  can  be  said  that  modulated  heating  and  cooling  of  the  airfoil 
surface  does  not  delay  the  formation  of  the  primary  recirculating  region  by  more  than 
2-3  degrees.  It  is  not  a  very  efficient  way  of  controlling  dynamic  stall.  Further  studies 
need  to  be  done  to  test  other  flow  control  techniques  like  generation  of  acoustic  pulses 
on  the  airfoil  surface  to  delay  the  formation  of  the  primary  recirculating  region. 
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surface  pressure 


Figure  79:  Surface  pressure  coefficient  at  a  =  18° 


Table  4:  CPU  and  Memory  Requirements  for  the  computations  on  Cray  C90 


Case 

Memory  (in  Mw) 

CPU-time  (in  hours) 

Ic 

9.8 

7 

7b 

128 

125 

9.14  Code  Performance 

The  Navier-Stokes  solver  is  fffily-vectorized  and  was  run  on  the  Cray  C90.  Some 
typical  CPU-time  and  memory  requirements  for  the  computations  are  shown  in  Table  4. 
The  solver  has  a  sustained  performance  of  approximately  380  MFLOPS. 
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10  CONCLUSIONS  AND  SCOPE  OF  FUTURE 


WORK 


10.1  Summary 

The  effects  of  compressibility,  pitch  rate  and  Reynolds  number  on  the  initial  stages  of 
boundary-layer  separation  on  a  NACA-0012  airfoil  pitching  about  its  quarter  chord  posi¬ 
tion  have  been  studied  numerically.  Computations  have  been  performed  using  two  sepa¬ 
rate  algorithms  for  the  compressible  laminar  Navier-Stokes  equations.  The  first  method, 
denoted  the  structured  grid  algorithm,  utilizes  a  structured,  boundary-fitted  C-grid  and 
employs  the  implicit  approximate-factorization  algorithm  of  Beam  and  Warming.  The 
solver  was  developed  originally  by  Visbal  [48,  49]  for  an  0-grid  topology  and  has  been 
modified  to  employ  a  C-grid,  Method  of  Characteristics  boundary  condition  at  the  outer 
boundary  and  the  Geometric  Conservation  Law  to  eliminate  the  grid  movement  related 
errors.  The  second  method,  denoted  the  unstructured  grid  algorithm,  utilizes  an  un¬ 
structured  grid  of  triangles  and  employs  the  flux-difference  splitting  method  of  Roe  and 
a  discrete  representation  of  Gauss’  Theorem  for  the  inviscid  and  viscous  terms,  respec¬ 
tively.  A  number  of  steady  and  unsteady  computations  have  been  performed  to  prove 
the  accuracy  of  the  two  algorithms,  including  the  boundary  layer  over  a  flat  plate,  sta¬ 
tionary  NACA-0012  airfoil  at  a  =  0°,  and  pitching  NACA-0015  airfoil.  The  test  cases 
were  compared  with  available  analytical  solutions  or  previous  computations.  The  effects 
of  implicit  and  explicit  artificial  dissipation,  and  the  temporal  and  spatial  resolution 
have  been  studied  in  detail  for  the  research  problem.  Accurate  grid-converged  solutions 
have  been  obtained  for  the  flow. 

The  principal  results  of  the  study  are  : 

•  The  emergence  of  the  primary  recirculating  region  has  been  traced  to  a  pair  of 
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critical  points  (a  center  and  a  saddle)  that  form  as  a  pure  shear  at  a  single  point 
and  immediately  bifurcate  and  move  apart. 

•  Secondary  and  tertiary  recirculating  regions  form  after  the  appearance  of  the  pri¬ 
mary  recirculating  region.  The  secondary  recirculating  region  interacts  with  the 
primary  recircidating  region  to  eject  the  fluid  close  to  the  wall  in  a  direction  ap¬ 
proximately  normal  to  the  wall.  The  ejection  of  the  fluid  near  the  wall  signifies 
boundary  layer  separation. 

•  The  accuracy  of  the  computations  has  been  confirmed  by  the  close  agreement  be¬ 
tween  the  structured  and  unstructured  grid  computations,  and  the  grid  refinement 
study  for  the  structured  grid  computations. 

Mach  number 

•  Increase  in  Mach  number  from  0.2  to  0.5  (at  fixed  Reynolds  number  and  pitch 
rate)  delays  the  formation  of  the  primary  recirculating  region  and  causes  it  to  form 
farther  from  the  airfoil  surface.  The  compressibility  effects  on  the  boundary  layer 
result  in  a  decrease  in  the  magnitude  of  the  leading  edge  suction  pressure  coefficient 
and  a  consequent  decrease  in  the  adverse  pressure  gradient  on  the  upper  surface. 
The  decrease  in  adverse  pressure  gradient  retards  the  movement  of  reversed  flow 
region  towards  the  leading  edge  and  hence  delays  the  formation  of  the  primary 
recirculating  region. 

Pitch  rate 

•  Increase  in  pitch  rate  from  0.05  through  0.2  (at  fixed  Mach  number  and  Reynolds 
number)  causes  the  transition  from  trailing  edge  separation  to  leading  edge  sepa¬ 
ration. 

•  Increase  in  pitch  rate  delays  the  formation  of  the  primary  recirculating  region  to 
higher  angles  of  attack.  This  is  due  to  the  decrease  in  the  adverse  pressure  gradient 
on  the  airfoil  upper  surface  which  retards  the  development  of  the  reversed  flow 
region  on  the  upper  surface  and  the  formation  of  the  primary  recirculating  region. 
The  primary  recirculating  region  also  appears  closer  to  the  leading  edge. 
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Reynolds  number 


•  Increasing  the  Reynolds  number  from  10^  to  10®  (at  fixed  Mach  number  and  pitch 
rate)  hastens  the  appearance  of  the  primary  recirculating  region  and  decreases  the 
length  scale  of  the  fiow  structures. 

•  Multiple  recirculating  regions  appear  simultaneously  near  the  leading  edge. 
Shock  wave 

•  A  shock  appears  on  the  top  surface  at  Rec  =  10®,  Moo  =  0.5,  and  12+  =  0.2  at 
a  =  19.5°.  A  separate  inviscid  computation  also  showed  a  similar  appearance  of 
shock,  implying  the  formation  of  the  shock  is  an  inviscid  phenomenon.  For  the 
inviscid  case,  the  shock  appears  closer  to  the  leading  edge  and  also  at  an  earlier 
angle  of  attack  as  compared  to  the  viscous  case. 

•  The  formation  of  the  recirculating  regions  near  the  leading  edge  was  found  not  to 
be  induced  by  the  shock. 

Stability  analysis 

•  Linear  stability  analysis  of  the  computed  flowfield  shows  that  the  formation  of  the 
primary  recirculating  region  is  related  to  the  convective  instability  of  the  flow.  It 
also  confirms  the  trend  seen  in  the  the  numerical  simulations  with  regards  to  the 
effect  of  change  in  pitch  rate  on  the  appearance  of  the  primary  recirculating  region. 

10,2  Future  Work 

The  field  of  unsteady  aerodynamics  is  relatively  in  the  nascent  stage.  A  substantial 
amount  of  research  has  recently  focused  on  the  study  of  the  complex  unsteady  effects  in 
the  fiow  over  a  pitching  and  oscillating  airfoil  so  that  ultimately  the  unsteady  effects  of  a 
pitching  airfoil  may  be  utilized  to  an  advantage  in  the  design  of  aerodynamic  bodies.  The 
present  work  presented  in  this  dissertation  has  concentrated  on  the  study  of  incipient 
stages  of  unsteady  boundary  layer  separation  and  the  effects  of  compressibility,  pitch 
rate  and  Reynolds  number  on  the  separation  process.  Computations  are  still  limited  to 
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relatively  low  Reynolds  number  flows  (with  respect  to  practical  range)  due  to  the  memory 
restrictions  of  the  computer.  Their  remains  scope  for  future  study  in  the  following  areas  : 

•  Practical  aerodynamic  flows  have  Reynolds  number  based  on  the  chord  of  the 
order  of  10®.  Increasing  Re,,  leads  to  more  expensive  computations,  both  in  terms 
of  memory  and  CPU.  Robust  and  efficient  algorithms  have  to  be  developed  to  solve 
such  flows.  The  use  of  parallel  computers  may  be  the  answer  to  the  problem. 

•  Flow  at  high  Re^  is  mainly  turbulent.  However,  there  are  not  any  turbulence 
models  which  can  accurately  simulate  the  flow  in  the  separated  regions  [6].  Proper 
turbulence  models  need  to  be  developed  for  such  flows. 

•  Study  of  3-D  effects  can  also  be  described  as  important.  But,  the  computer  re¬ 
sources  available  at  the  present  time  does  not  allow  an  accurate  simulation  for  3-D 
flowfields. 

•  The  control  of  dynamic  stall  process  is  a  very  important  problem  that  needs  to 
be  addressed.  In  this  study,  the  formation  of  the  dynamic  stall  vortex  has  been 
traced  to  a  pair  of  critical  points  in  the  flow  field.  It  might  be  of  interest  to  apply 
different  flow  control  techniques  at  the  location  of  the  formation  of  these  points  to 
ultimately  delay  the  formation  of  the  dynamic  stall  vortex. 
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A  Jacobian  Matrices  for  Structured  Grid 

Algorithm 

Explicit  form  of  the  Jacobian  matrices  has  been  presented  in  this  appendix.  The 
matrices  are  obtained  analytically  by  rewriting  the  matrices  F,  G,  Vi,  V2,  and  W2 
in  terms  of  p/J,  pu/J,  P'^/ J-,  and  pej  J  and  differentiating  according  to  eqn.  (21)  and 
eqns.  (25)  to  (29). 

The  Jacobian  matrices  A  and  B  given  by  eqns.  (21)  and  (25),  respectively  are  : 

^6  ix  4  0  ^ 

ix<i>-uU  U-{'i-2%u  £y«-(7-l)^®v  (7-1)^® 

^y(f>-vU  ^xV  -  {7  -  lUyU  U-{'i-2)iyV  (7-1)^!/ 

^U{24>-'ie)  {'ie- -l)uU  (76  -  ^)^y  -  (7  -  l)vCf  7!^ +  6  y 

where 

U  =  U-U-,  V  =  v-Vf,  't’=\('t-^^)(v?+v^)  (109) 

U  and  V  are  the  contravariant  velocities  described  by  eqn.  (7). 

Tjt  Tjx  T)y  0 

T]x(l>  -uV  V  -  (7  -  2)t]xU  TtyU  -  (7  -  l)»7®t'  (7  -  1)^® 

rjyCj)  -  vV  T]xV  -  (7  -  l)»?y«  V  -  (7  -  2)T]yV  (7  -  l)T]y 

V{24>  -  7e)  (7e  -  ^)j7x  -  (7  -  1)«1^  (7^  -  <t>)Vy  -  (7  “  '^)vV  -yV  +  Tjt 
Equations  (27)  and  (29)  define  the  Jacobian  matrices  R  and  S,  respectively. 

0  0  0  0 

{biu  +  b2v)  61  62  0 

{b2U  +  63^)  62  63  0 

Rjll  Ra2  Riz  RiA 
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where 


Rii  =  —{hiU^  +  2h2uv  +  +  647(7  “  +  v^  —  e)} 

Ra2  =  {61  -  647(7  -  1)}“  +  ^2V 
R43  =  h2U  +  {63  -  647(7  -  1)}^’ 

R44  =  647(7  -  1) 


where 


S 


0  0 

1  -{dxu  +  d2v)  di 

P  —{d2U  +  dzv)  dz 

L  iS'41  S42 


0  0  ^ 
<^2  0 
dz  0 

^43  544  ^ 


(112) 

(113) 

(114) 

(115) 


(116) 


‘S'41  =  —{dxu^  +  2d2uv  +  dzv^  +  <^47(7  —  l)('w^  +  f  ^  —  e)}  (117) 

S42  =  {<^1  -  <^47(7  -  l)}w  +  dzv  (118) 

Siz  =  dzu  +  {dz  -  <^47(7  -  l)}'i^  (119) 

S44  =  d4'y{'y  -  1)  (120) 


Jacobian  matrices  P  and  R  are  defined  by  eqns.  (26)  and  (28).  The  matrices  (— P  + 
R^)  and  {—Q  +  Sr,)  required  in  the  Beam- Warming  algorithm  (eqns.  (39)  and  (40))  are  : 
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coefficients  hi  (i=l,....4)  are  given  in  eqn.  (8). 


coefficients  di  (i=l,.,..4)  are  given  in  eq.  (10). 


B  Geometric  Conservation  Law 


Geometric  Conservation  Law  is  based  on  the  basic  principle  that  a  uniform  steady 
flow  in  a  computational  domain  should  remain  uniform  if  there  are  no  external  distur¬ 
bances  and  even  if  the  computational  grid  moves  with  time. 

The  unsteady,  compressible,  two-dimensional  Navier-Stokes  equations  in  strong  con¬ 
servation  form  has  been  presented  in  eqn.  1.  The  corresponding  equations  in  transformed 
coordinates  (^,  Tj)  can  be  written  as  : 


and 


^a:  —  VriJ]  "b 

Vx  =  -y^J ;  Vy  =  ;  Vt  =  {xty^  -  ytx^)J 

Beam- Warming  algorithm  is  applied  to  eqn.  131  with  trapezoidal  time  differencing  to 
give  ; 

Ag"  =  Y  +  ‘5)]”(137) 

where 

Ag"  =  -  g"  =  A  0)  =  AU”  -k  U^A  0)"  (138) 

no 


Substituting  eqn.  138  in  eqn.  137  gives  : 


-A<[|(^-K)  +  A(e_5)]”-c;»A(i)” 

Approximate  factorization  of  eqn.  139  leads  to  : 


At  d 


2  dt 


At  d 


2  a??' 


where  B,  C,  and  V  are  the  Jacobian  matrices  defined  by  : 

B  =  c  =  ^-  P=^ 

dq'  dq’  dq'  dq 


(139) 


(141) 


If  we  consider  uniform  steady  flow  in  the  computational  domain  then  the  flow  should 
remain  unaltered  after  a  time  At.  Also 


K,j  =  «oo;  vfj  =  Voc]  vli  =  Poo;  plj  =  poo;  u  =  Constant 
A{7”  =  0;  F  =  Constant]  G  =  Constant]  R  =  [0];  S  =  [0]  (142) 


Substituting  eqn.  142  and  eqns.  133  to  136  in  eqn.  140  leads  to  : 


The  metrics  a:^,  and  yr,  are  calculated  numerically  using  finite-difference  formulas. 
Figure  80  shows  9  nodes  of  a  part  of  the  computational  transformed  grid.  If  central  finite- 
differences  are  used  to  compute  the  metrics  at  point  ‘C’  of  the  grid,  and  A|  and  Arj  are 
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assumed  to  be  1.0,  then  : 


dyr,  dy^ 

d(  dy 


dx^  dXr, 
dr]  di 


""  2^^^^  "  ^  (5)  “  ^1  “  ^8  +  y&)  - 

(^)  (2/3  -  2/8  -  2/1  +  2/6)  =  0 

-  2:4)  -  -  *7)  =  Q)  (2:3  -  ®8  -  a:i  +  ®6)  - 

(^)  -  *1  -  *8  +  aJe)  =  0 


(146) 

(147) 


Figure  80:  Part  of  the  transformed  (^,  y)  grid. 
Substitution  of  eqns.  (146)  and  (147)  in  eqn.  (145)  gives  : 


For  a  stationary  grid  A(l/J)  is  equal  to  0.  Geometric  Conservation  Law  states  that 
for  a  grid  moving  with  time,  the  term  —U'^A{lfJ)'^  (given  by  eqn.  (150))  is  non-zero 
and  has  to  be  appended  to  the  right  hand  side  of  the  Beam- Warming  algorithm. 
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C  Method  of  Characteristics 


The  Method  of  Characteristics  is  based  on  the  Euler  equations.  Figure  81  shows 
the  outward  moving  perturbation  wave  generated  by  a  point  source  after  a  time  At. 
There  is  a  uniform  flow  of  velocity  C/qo  coming  into  the  domain.  The  perturbation  wave 
moves  in  an  approximately  elliptical  pattern.  If  the  outer  boundary  is  of  the  shape  of 
the  perturbation  wave,  the  solution  can  be  considered  to  be  one- dimensional  near  the 
outer  boundary.  A  C-grid  matches  with  the  wave  shape  to  a  certain  degree,  when  the 
outer  boundary  is  sufficiently  far  away  from  the  airfoil.  Therefore,  a  one- dimensional 
Method  of  Characteristics  can  be  used  at  the  outer  boundary  of  the  C-grid.  In  the 
Method  of  Characteristics,  the  Riemann  variables  corresponding  to  the  characteristics 
moving  out  of  the  domain  is  interpolated  from  the  inside  of  the  domain  and  all  the 
Riemann  variables  corresponding  to  the  characteristics  moving  inside  the  domain  are 
set  as  boundary  conditions  at  the  computational  outer  boundary. 

A  local  coordinate  system  {x,  y)  is  constructed  orthogonal  to  the  outer  boundary 
with  X  normal  to  and  directed  outward  from  the  boundary.  The  derivatives  along  the 
boundary  (y)  are  neglected  (one-dimensional  problem).  Therefore,  )  =  0. 

The  continuity  equation  can  be  written  as  : 


Dp  du 


(161) 


where  u  is  the  velocity  in  x-direction.  Also  p  =  p{p,  s),  where  p  is  the  pressure  and  s  is 
the  entropy.  Therefore, 


r  f ^P\  J  f ^P\  J 

\  ^  /  3  \  /  p 


(162) 


For  locally  isentropic  flow,  ds  =  0.  From  eqns.  (151)  and  (152)  : 

Dp  / Dp  1  Dp  _  du 
Dt  \dp) ^  Dt  a?  Dt  ^ dx 


(153) 
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outer  boundary 


perturbation  wave 


Jairfoil  (source) 


(a-U)At  (a  +  U)lt 


Figure  81:  Perturbation  wave  around  an  airfoil. 
1  /  ^9p\  du 

The  momentum  equation  can  be  written  as  : 

du  _du  1  dp 
dt  '^^dx^  pdx~ 

Addition  of  eqns.  (154)  and  (155)  gives  : 

du  5ul  ^  \ dp 

=0 

Subtraction  of  eqn.  (154)  from  eqn.  (155)  gives  : 

/_  ^  \dp  -.dp] 

^  +  (»-a)g=J--[-  +  («-a)-J=0 

Let,  ^  —  u  +  a.  Therefore, 

du  du^  _du  du 

dt  dt  dx  dt  dt  ^  dx 


dp  dp  dp  dx  dp  ,  .  dp 

dt  dt  ^  dx  dt  ^  ^  dx 


(160) 


Substitution  of  eqns.  (158)  and  (159)  in  eqn.  (156)  gives  : 


For  a  perfect  gas 


+  —  =  0 


pa 


a^  =  — 


(161) 


Also,  for  an  isentropic  process 


p  = 


(162) 


where  C3,  and  C4  are  constants. 


dp  =  C4 


r-1/ 


(163) 


and,  p  =  C4')a^^ 


(164) 


From  eqns.  (160),  (163),  and  (164) 


_  f  dp 

u+  —  =  constant 
J  pa 


(165) 


or,  u  + 


2a 

7-1 


=  constant 


Therefore,  the  first  characteristic  equation  is 

2a 


dx 


u  H - -  =  constant  along  —  =  u  +  a 

7  —  1  dt 


Now,  let  ^  =  u  —  a.  This  gives  : 


du  du  dudx  du  du 

dt  dt  ^  dx  dt  dt  ^  ^  dx 


(166) 


(167) 


(168) 


dp  dp  dp  dx  dp  dp 

dt  dt  dx  dt  dt  dx 

Substitution  of  eqns.  (168)  and  (169)  in  eqn.  (157)  gives  : 

*-^  =  0 

pa 


(169) 


(170) 
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(171) 


—  =  constant 
pa 


u  — 


2a 

7-1 


=  constant 


The  second  characteristic  equation  is 


_  2a  ,  dx  _ 

u - =  constant  along  —  —  u  —  a 

7  —  1  dt 

The  energy  equation  can  be  written  as  : 


Ds  ds  _ds 

Let  ^  =  u.  Then 

ds  ds  ds  dx  ds  _ds 
dt  dt  ^  dx  dt  dt  dx 

From  eqns.  (174)  and  (175) 


The  third  characteristic  equation  is 


s  =  constant 


along 


dx 

dt 


=  u 


(172) 

(173) 

(174) 

(175) 

(176) 

(177) 


The  equations  are  assumed  to  be  one-dimensional  at  the  boundary.  Therefore,  the 
velocity  in  the  y-direction  (tJ)  can  be  assumed  to  be  constant  along  the  characteristic 

dt  “• 

The  fourth  characteristic  equation  is 


V  =  constant 


dx  _ 
along  —  =  u 


(178) 
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D  Jacobian  Matrices  for  Unstructured  Grid 

Algorithm 


D.l  Inviscid  Jacobian 

Consider  the  terms  in  the  inviscid  Jacobian  (90)  involving  the  matrices 

2^  dQl 

A  =  -T~'^^^As 
"  "  dQr 

The  matrices  Ai  and  Ar  have  an  identical  form,  differing  only  in  their  evaluation  using 
Qi  and  Qr,  respectively.  Define  the  matrix  A{Q)  by 

A  =  T-'^^As 


Note  that  A  does  not  include  the  leading  factor  of 
The  components  Aij  are 


Ai2 

Ai3 

Ai4 

A21 

A22 

A23 

A24 

A31 

A32 


—UgAs 

Ay 

—Ax 

0 

2  [(7—3)  +  (7  — 1)  Ay  +  uvAx 

(3—7)  uAy  —  vAx  —  UgAs 
—  (7  —  1)  vAy  —  uAx 

(7-1)  Ay 

— uuAy  +  I  [^(3— 7)u^  —  (7— l)u^  Ax 
vAy  +  (7— 1)  uAx 
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=  uAy  —  {3—^)vAx  —  UgAs 
A34  =  -(7-1)Ax 

A41  =  — 7146  +  (7  — 1)  Ay  +  7t;e  —  (7  — 1)  Ax 

A42  =  ~  "^”2 — (3'*^^  + 'y^)|  Ay  +  (7— 1)  wt)Ax 

A43  =  —  (7  — 1)  ttwAy  —  76  —  ^  ^  ^  +  3i;^j  Ax 

A44  =  7  {uAy  —  vAx)  —  UgAs 


where 


Ay 

Ug  =  Ug- - Vg 


Ax 

As 


Next,  consider  the  terms  in  the  in  viscid  Jacobian  (90)  involving  the  matrices 


Bi 

Br 


where 


C  =  lT-^S\A\S-'^TAs 

Define  the  vector 

Q  =  (w,v,  a)^ 

Note  that  Q  has  three  components.  Then 

dC  _  dC  dQ 
dQi  dQ  dQi 
dC  _  dC  dQ 
dQr  dQ  dQr 

where  dC / dQ  is  a  third  order  tensor 

^  =  \  T-'^DTAs 
dQ  ^ 

and 

£(£|^ 

dQ 
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The  components  Dijk  are  defined  as 

a(s|A|s-'|y) 

^iok  -  - Ta - 

dQk 

where  |tj  indicates  the  component  in  the  row  and  column  of  the  matrix 
The  expressions  are 

-Dm  =  ^ l^sl  +  IA4I  —  2IA2I  +  As  +  A4  —  2A2  +—  IA4I  —  l-^sl 

2  a‘^  ^  2  ^  2a  '■ 

U  \r  r  1  r 

+  A4  —  As  +  A2 

2a 

D112  =  |As|  +  IA4I  —  2IA2I 

Diis  =  —(7  —  1):^  |As|  +  IA4I  —  2IA2!]  +  As  —  A4  — IA4I  —  [Asl 

CL  ^  CL  ^  CL 

-1|1a..A3] 

Di21  =  ^-2  [^1^2!  —  IA3I  —  IA4I  +  2  ^  ^^2  —  As  —  A4  As  —  A4 

D122  =  0 

Di23  =  -  (7-I)  ^  2IA2I  -  |As|  -  IA4I]  +  A4  -  As  -  —  JAsI  -  IA4I 

1  ri  I  1 

+  rT  AS  +  A4 

2iQ/ 

'i  —  Ivr*  *  ~  ^ 

Disi  =  —1  [2A2  —  As  —  A4J 

Zt  CL 

Di32  =  [2IA2I  —  |As|  —  IA4II 

Dis3  =  ~  (T"!)  ^  [2IA2I  -  IA3I  -  IA4I]  +  A4  -  As 

Di4i  =  As  +  A4  —  2A2 

Di42  =  0 

Di4s  =  —  ^3-  |As|  +  IA4I  —  2IA2I  +  [-^3  —  A4 

D211  =  |As|-^— :r  +  2^0  +  -  |As|;;^  [2^  +  a]  +  |A4|^^^— ^  -  2'ua  + 

4a^  2a 

+  IA4I— [2«  —  a]  +  IA2I  1  —  ^_2  +  2u  j  +  A2  —  ^^^2 . ^ uq  +  u 


g(a  +  a)  -  —  (S  +  a)  -f  A4  -  a)  +  —  -  a) 


D212  —  [lAsKu  +  a)  +  |A4|('a  —  a)  —  2|A2|a 

Z  CL 
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D 


213 


D222 

D223 


D231 

D232 

D233 


D24I 

D242 

D243 

Dsw 

D312 

D313 

D321 

D322 

D323 


il  _ 

2a?  2a? 


{a  +  2u)  q 


+  IA4I 


+  IA2I 


(7~1):i/-2  ,  a2^ 

— — — u{u  +V  ) 


a*' 


+  A3{u  +  a) 


(7  — 1), 


_i_  I  (lull  o 

2a  2a?  ^ 


-A4{u  —  a) 


u  (7  — 1)  - 

+— +  M:rT^g 


D221  —  +IA2I 


2(7-1) 


2a  2a? 
u 


a 


+A2 

0 

-IA2I 


ft-4 


+  IA3I 

+  A3('U  +  a) 


1  (7— 1)  U 


-|A4i 


’r  +  47^(2fi-S) 


2a 


2  -2 


a 


2a  2a? 
A4{u  —  a) 


1  ,  (7-I)  u 
2a  2  5? 


2(t-1)| 


+  IA3I 


~2  ~ 

(7-i)|pr  +  ^  “i? 


+  A3-<®  +  “' 


2a 


1  -  (7-1) 


+  IA4I 
(7-1)  ^ 


^2  ; 

(7-1)  I  +^2 


+  A, 


(it  —  a) 


2a 


u 


1  +  (7-1)t 


—  [2IA2I  +  2'uA2  —  IA3I  —  (it  +  a)A3  —  IA4I  —  (u  —  a)A4 

a 


y--^j 

2a2 

2^|A2|  -  {u  +  a)|A3|  -  («  -  a)|A4l 

2IA2I 

+  IA3I 

7“1  ^  /«-  -X 
.3(2«  +  a) 

^  d 

+ 

1  1 

’7-1  ^  //,=  ~x' 
2 

—A3 

7-1 

2a? 

0 

-IA3I 

+A3 


7-1  V 

^(«  +  a) 

^  CL 


+  A4 


7-1^/^  -A 


IA3I  +  IA4I  —  2IA2I  —  2uA2{u  +  a)A3  +  (iZ  —  a)A4j 

+  iA2| 


1^(26  + a) 


7-1/^  I 

-^(x  +  a) 


+  1^4 

.  r 
-A4 


2(7-1)^ 


7-1/=  ~x 
— r(«  -  «) 


2a" 


t;(I>iii  -  Ai) 

(7-1) 


2a2 

^^-^113 


IA3I  +  IA4I  —  2IA2I  ^  IA4I  —  lAslj  +  vH 


'112 


vDi21 

7  — 1  u 
~2 
vD\23 


ui  r  /V  ^  -^1 

T2  2IA2I  —  IA3I  —  IA4I  +—  IA3I  —  IA4I  + 'ijl?l22 

O'-  ^  za  '■ 
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s*  I  Si< 


Ai  +  vD\3\ 

"^"2  —  I'^sl  —  IA4I  -\-vDiz2 

vD\3z 


vDxil 

^~2  1-^31  +  1^41—213(21  +  ^-0142 

V-D143 


fi|(,-.)i-s|} 


it  I  {(H  —  ua)  \.  .u  1  (tt  — a)  r,  .q 

+i'+Sr^  (^-'>1+1  +Sr^+-i)!+“ 


+  |A2|«  1-2(7-1)4 

Qj  ,  2i€L 

r  (H  —  Ua)  r,  .  T  r  /  V 

+A4^^ — TT — -  (7-l)T  +  ti  +A2  -{'y-1)^  +  q-v‘ 

JdCL  .CL  ,  CL 

(7  — 1)  V  r  ~  i~x|T  I  o~|r  I  ,  irr 


(7-l)r-« 
L  a 


rj  [{H  +  Sa)|A3|  -  2?|A2|  +  {H  -  M)\U 

CL 


+^f  [l^3l  +  |A,|  -  2|Aj|]  +  II  [|A,|  -  IA3I]  -5IA2I 

+IA3I  (+4^  f-(7-l)l  +  |]  +  4  [^2^1  -  I 

[  a  2j  \7— 1  /  ®  L  ^  a  2 

+IA,|  1+^  _(^_i)|  _  !  +  f  J!_  _  s')  i  flizi) i  +  ! 

\  [  a  2j  \7  — 1  )  a  [  2  a  2 

f  (F  +  «a)  r(7-l)g  Sl'l  ~  f(#-Sa)  [(7-1)9  ,  -s' 

a  ^a”2  r  ‘1  i  2“a  +  2 


9  _  «]  1  _ 

'd  2J  J 


+  l'^2| 

-1+ 1 (1^+^  ,  [  ,  (,-i)i] }  ,  |A3|(i^  (,+ 

~  (^  +  tZa)  [^  ,  ^.u]  ~  (H  —  ud)  [^  .  ^.u]  ~  .  _^.uq 

2a  P  -  ("-"aj  -  2H  P  +  (^-')aj  + 

[2|A3|-|A3|-|A4|]+^[|A,|-|A4|] 

.  {(H  +  ud)  \.  _iZ  l1  [2a  ~1  1  [,  ,  ..mII 

+W  (^-i)a-2  +  —1+^  ra  '-("-I'a  | 
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Dizi 

D432 

D433 


D441 

-D442 

D443 


where 


+  IA4I 


—  ua) 


+A3 


{H  +  ua) 


/  1 
(7-1)t +  - 
a  I 


2a 


—  u 


2d 


u 


a 


+  A4 


7-1 

{H  —  ua) 
2d 


2d 


u 


1  +  (7-I)- 


1  +  (7-I) 


u 


“IA2I 

(7-I)  V 

2  d? 

{H  —  ud)A.4  +  vA 


“IA3I 
+  IA2I 


(7-1) 

2d^ 

(7-1) 


[2tt|A2|  +  25A2  -  (tt  +  «)|A3|  -  {H  +  ud)A3  -  {it  -  a)|A4|- 

(7-1) 


{H  +  ud)  +  v  —  IA4I 


2d^ 


[H  —  ud)  + 


+  +  1 


a 


'aM 

(H  —  ud)  1 


2d 


2  [7-1 


+  u 


+  |A4|(7-1): 


2d 


—As 


d  2  [7— 1 
(7-1)  V 


—  u 


-21A2|(7-1)P 


2  dr 


H  -\-ud 


}  +  A,{ 


(7-I)  ^ 


H  —  ud 


2  d^ 

[(S  +  a)|A3|  +  {H  +  Sa)A3  +  (S  -  a)|A4|  +  {H  -  iid)A4  -  2S|A2|  -  2qA2 

7— 1  V 


2  d 


3I  4-  IA4I  —  2IA2I 


+  IA3I 


(7-1)  /  {H  +  iid)  1 


{H  —  ud)  1 
d  ^  2 


a 

2d 


+  - 


2d 


2  [7-1 


+  u 


+  |A4|^ 


[7-1 


—  u 


+  |A.I(7-l)f  +  A3i^fc®“) 


a=^ 


~  (7-l)(F-tta) 
-A4 — - — 


q  =  \{u  +v^) 

The  terms  Aj  =  ||  Aj  |1  jdXj  are 

^  _  sign(Aj)  if|Aj|>5 

"  \Xj/S  if|A,l<5 

The  scalar  quantities  Aj,  j  =  1,2, 3, 4,  should  not  be  confused  with  the  diagonal  matrix 
A  in  (75). 
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£5.  = 

w  -  s®r3  - 

EU  =  W  -  sb;,  -  S£;j 


where 


(Tl  = 

2^  PI 

-  [K 

Gr  = 

2  y  Pr 

Q  = 

K[  = 

1  f(7-l) 

p/  +  ^  \  2 

K[  = 

(7-l)«i 
{pi  +  q) 

Ki  = 

{'y-l)vi 
{pi  +  q) 

K[  = 

7 

{pi  +  q) 

k:  = 

1  f(7-l) 

1 

Pr  +  ^  \  2 

Kl  = 

(7-l)«r 

{Pr  +  Q) 
{l-l)Vr 

k;  = 

{pr  +  Q) 

Kl  = 

7 

{pr  +  Q) 

2  ,  2 
■**/  +  Vj 


^  +  i?/ 


2  ,  2 
«r  +^r 


The  inviscid  Jacobian  is  completed  by  specification  of  the  matrices  dQijdQj  and 
dQrfdQj  in  (90).  These  matrices  follow  directly  from  (77)  and  (78).  Assuming  linear 
reconstruction  of  the  conservative  variables  (Fig.  9), 


dQi 

dQj 


I  if  j  =  i 

<;a^Q  if  j  is  point  a 
if  j  is  point  b 
if  j  is  point  c 


124 


where  I  is  the  identity  matrix  and 


and 


ft 


ft 


■z^{nb  +  fie)  •  fds 
^Kbc 

— — (tcc  +  fia)  •  rds 
^Vabc 

(n„  +  fih)  •  fds 

^Vabc 


^p«?  ^pvt  ^pe) 


The  terms  $p,  $pu,  and  $pe  represent  the  limiter  function  $  for  the  conservation  of 
mass,  X—  and  y—  momentum,  and  energy,  respectively. 

A  similar  expression  for  dQrjdQj  may  be  obtained  taking  into  account  the  location 
of  the  cells  employed  to  form  Q^. 

Special  modificatioiis  are  required  for  boundary  cells,  and  are  not  presented  herein. 


D.2  Viscous  Jacobian 

The  viscous  Jacobian  (92)  may  be  expressed  as 

dN  _  dN  dQd 

dQj  ~  dQadQj^  dQbdQj^  dQcdQi^  dQd^Qj 

where 

Q  = 

and  the  subscripts  a  tod  refer  to  the  points  shown  in  Fig.  3.  Define 

dN 

dQa 

and  similarly  for  andTV'^.  Then 

=  0 

^2  =  0 
A“3  =  0 

1^2“  =  Eoo(-|A7/A|/4i  -  Aa:Aa;4i) 

N22  =  Hoo(-|AyAa;4i  +  Aa:A?/4i) 
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7V“ 

-‘''23 


=  0 


?V“ 

■‘''31 


■'''41 

-'''42 

■‘''43 


-‘''ll 


E:oo(Ai/As4i  -  |Aa;Ai/4i) 
-E:oo(AyAy4i  +  |Aa:Aa:4i) 

0 

+  Ud)N^^  +  \{vb  +  Vd)N^^ 

^{Ub  +  Ud)N22  -f  \{Vb  +  Vd)Nt2 

■(Aj/Ai/4i  +  Aa:Aa;4i) 


Pr(7  — 1) 


AT**  — 
-‘''12  — 


TV**  — 

-‘''13  ~ 


0 

0 

0 


•A^oi  =  -Soo(|AyAi/i2  +  Aa:Axi2) 
Soo(-|Ai/Axi2  +  AxAl/12) 


'21 
NI2 
Nl, 

Mb 

■‘''31 

■‘''32 


33 


= 

2/1  ' 

=  Soo(Ai/Axi2  -  |AxAyi2) 

=  -EociAyAyu  +  |AxAxi2) 

=  -^N, 

2ndT  ^ 


■^41  ~  2  ^d)^21  +  2('^b  +  Vd)N2i  +  ^N2 

^42  ~  2('^b  +  y'd)N22  +  \{Vb  +  Vd)N^2  2^^ 


#  _ 

■‘''43  ~  -  — 


N4  -  lN2{Ub  +  Ud)  -  lN3{Vb  +  Vd)  - 


2/1  dT 

+  |(^ib  +  l^ci) -^23  +  2(^‘’ 


(AyAyi 


■‘’11 


Arc 

■‘V12 


A/c 

■‘'^13 


0 

0 

0 


Arc 

■‘''21 


-Hoo(|Ai/A?/23  +  AxAx23) 


^22 


,(-|At/Ax23  +  AxAt/23) 


^23 


=  0 


Arc 

■‘V31 


Eoo(A?/Ax23  -  \AxAy2z) 


+  AxAxi2) 
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^32 

= 

NI3 

= 

= 

-‘*42 

= 

■‘''43 

= 

Nil 

= 

-‘''12 

= 

-‘’13 

= 

-‘’21 

= 

= 

-‘’23 

= 

= 

^32 

= 

-‘’33 

= 

-‘’41 

= 

-‘’42 

z=: 

-‘’43 

= 

where 

-E:oo(Ai/A?/23  +  |Aa:Aa:23) 

0 

l{Ub  +  Ud)N21  +  1(^6  + 

+  Ud)N22  +  \{Vb  +  Vd)N^2 
1 

-{AyAy23  +  AsAaija)  +  +  'Wd)A^23  + 


Pr(7  — 1) 


0 


0 
0 

-Soo(|At/At/34  +  AxAx34) 

Eoo(-|At/A®34  +  A®Ay34) 

2iJ.dT  ^ 

=  Eoo(AyAx34  -  |AxA2/34) 

=  -Soo(A2/Aj/34  +  |AxAx34) 

2/i(iT  ^ 

+  Ud)N2i  +  \{Vb  +  Vd)N2i  +  |iV2 
\{ub  +  Ud)N22  +  +  Vd)N^2  "I"  2-^3 

N4  -  lN2{ub  +  Ud)  -  lN3{vb  +  Vd)] 

1 


1  dy, 
2fi  dT 


’  Pr{'y  —  1) 


{AyAy34  +  AXAX34)  +  |(«i>  +  Ud)N23  4 


M(X)y 


2ReooV^ 


and  Vu  is  the  volume  of  the  quadrilateral  defined  by  point  a,  b, 


Axi  =  Xb  —  Xa 
Ax2  =  Xc  —  Xb 
Ax3  =  Xd  —  Xc 
Ax4  “  Xd  Xd 

Ayi  =  yb-ya 


2  ('*^6  4  Vd)N^2 


livb  +  Vd)N22 


c,  d  in  Fig.  3  and 
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At/2  =  yc-yb 
At/3  =  yd-yc 
At/4  =  ya- yd 

while 


Aa:i2  =  AaJi  +  Ax2 
Aa:23  =  Ax2  +  Axs 
Ax34  —  Axs  +  Ax4 
Ax4i  =  A«4  +  Asi 


and  similarly  for  Ayi2,  ....  Note  that  Ax  and  Ay  refer  to  the  x—  and  y—  spacing  on 
edge  k. 

The  matrix  dQa(dQj  may  be  expanded  as 

dQa  _  dQa dQa 

dQj  dQa dQj 

V  A 

Define  Q  =  dQtfdQt  where  the  subscript  ,  implies  points  a,  b,  c  or  d.  Then, 


P 


Qi2  -  - 
P 

Qiz  =  0 


Qi4 

Q21 

Q22 

Q23 

Q24 

Qzi 

Q32 

Q33 

Q34 


0 

V 

p 

0 

1 


P 

0 

(-pe  +  pv?  + 

-7(7-1)- 

P 

-7(7-1)^ 

P 

7(7-1) 

P 
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where  the  terms  are  evaluated  at  the  respective  point.  The  remaining  matrix  is 


/  if  •  and  j  correspond 
to  the  same  centroid 
if  •  corresponds  to  a  node 

where  I  is  the  identity  matrix  and 


dQ, 

dQj 


^  ^Ifor  cell  j 

=  “F - - 

^2cells 

Special  modifications  are  required  for  boundary  cells,  and  are  not  presented  herein. 


D.3  Roe  Flux 

For  completeness,  the  elements  of  the  Roe  flux  term  S\A\S~^AR  in  (75)  are  provided. 


'^11 

= 

0 

5i2 

= 

1 

5i3 

= 

1 

Sl4 

= 

1 

S21 

0 

S22 

= 

u 

S23 

= 

u  +  a 

§24 

= 

u  —  d 

z= 

V 

§32 

= 

V 

§33 

z= 

V 

§34 

= 

V 

§41 

= 

~2 

V 

§42 

= 

q 

§43 

= 

H  +  ud 

§44 

= 

H  —  ud 
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|A|  = 


The  vector  AR  =  Ri  —  Rr  where 


-1 


(7-1)  ~ 


a 

(7-l)| 

7-1 


a 

(7-1)-  u 

2a?  ^  2a 
(7-1)  u  1 

2  2a 

(7-1)  ^ 


a^ 


7-1 

2a2 

(7-1) 

2a 


u 


(7  — 1)  U  1 

2  2a 

(7-1)  V 

2  a? 

7-1 

2a? 

diag(|Aa|,|A2|,|A3UA4|) 


R  =  {p,pu,pv,pef 
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E  Stability  Analysis 


This  appendix  formulates  the  equations  for  a  linear  hydrodynamic  stability  anal¬ 
ysis.  The  non-parallel  Orr-Sommerfeld  equation  has  been  derived  in  a  rotating  frame  of 
reference  assuming  incompressible  flow  and  constant  pitch-rate.  The  boundary  condi¬ 
tions  have  also  been  derived  for  a  constant  rate  pitching  airfoil. 

E.l  Governing  Equation  for  Stability  Analysis 


The  non-dimensional  continuity  and  momentum  equations  for  an  incompressible  flow 
in  a  reference  frame  attached  to  the  airfoil  pitching  at  a  constant  pitch  rate  Q,  can  be 
written  as 

ri'ii.  rill 

(179) 


du  P 

dx  dy 


du  du  du 


Xdp  1 

2^1;  — - — - h 

pdx 


dv  dv  dv  \dv  1 

+  u—  +  V—  -t-  -f 


(180) 


(181) 


dt  dx  dy  pdy  Re  \dx^  '  dy^^ 

A  small  disturbance  is  introduced  into  the  base  flow  an  the  development  of  the  flow 
field  in  time  is  monitored  to  determine  whether  the  flow  is  stable  or  unstable.  The 
velocities  and  pressure  are  written  as  follows 


u  =  U  +  u' 

(182) 

V  =  V  -\-v' 

(183) 

p  =  P  +p' 

(184) 

These  expressions  are  substituted  into  equations  (179), 

linearized  to  give  the  following  equations 

(180),  and  (181)  and  then 

du'  dv' 
dx  ^  dy 

(185) 
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du'  ,dU  ,dU  Bu'  ^Bu'  , 

— t"  U  —  h  V  — - 2fiv 

Bt  Bx  By  Bx  By 


1B£  J_  /B^u'  B\' 

p  Bx  ^  Re  \  Bx^  By^ 


Bv'  ,BV  ,BV  ,,Bv'  ^Bv'  ,  IBp'  \  ( B^v'  BH' 
■ft  +  “  + ’' aj;"  +  ^ V 

The  velocity  and  pressure  perturbations  may  be  expressed  as 


(186) 


(187) 


u'{x,y,t)  —  u{y)  exp{i{kx  —  a;t)} 


(188) 


v'{x,y,t)  =  v{y)  exp{i(A:x  —  a;t)} 
v'i^tVtt)  =  p(y)  exf{t{kx  -  Ojt)} 


(189) 

(190) 


where  k  and  u)  are  the  complex  wavenumber  and  frequency  respectively,  and  u{y),  v(y) 
and  p(y)  are  the  complex  complex  amplitude  functions  for  the  velocities  and  pressure. 
Substituting  equations  188,  189,  and  190  in  equations  185,  186,  and  187,  we  get 


iku  +  Dv  =  0 


(191) 


—uijji  +  +  Uiiki  +  Vu  —  2Clv  =  —pki  +  -^{—uk^  +  D^u)  (192) 

Bx  By  Re 

1. 

—vivi  +  u— — h  V— — h  Uvki  +  Vv  +  2Qu  =  —Dp  +  —(—vk^  +  D^u)  (193) 
Bx  By  Re 

where  “D”  denotes  differentiation  with  respect  to  y.  The  velocity  perturbations  can  be 
expressed  in  terms  of  complex  streamfunction  perturbation 


,  Bf 


(194) 


where 


=  viy)  exp{i();3:  -  u>i)} 


(195) 


(p{y)  is  the  complex  amplitude  function.  Utilization  of  equations  194  and  195  eliminates 
equation  191.  From  equations  192  and  193  we  get 

i-iojD  +  +  iUkD  +  VD^  -  i^k  +  2ifcO  )  ip  =  -ikp^^{-k'^D  +  D^)ip  (196) 

\  Bx  By  j  Re 

( -few  -  ik^  +  k^U  -  ikVD  +  +  2^0]  if  =  -Dp  +  -^{ik^  -  ikD^)ip  (197) 

\  By  Bx  J  Re 

From  equation  196 

.  fijD  .BUD  .VD^  .kD  .D^  BU 
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The  above  equation  is  the  Orr-Sommerfeld  equation  for  a  non-parallel  flow.  It  is  a 
fourth  order  complex  differential  equation  which  requires  2  boundary  conditions  at  each 
of  the  two  boundaries  to  solve  the  eigenvalue  problem.  The  following  subsection  presents 
the  appropriate  boundary  conditions  for  a  constant  rate  pitching  airfoil. 


E.2  Boundary  Condition  for  Stability  Analysis 

Suitable  boundary  conditions  need  to  be  specified  at  the  two  boundaries  to  solve  the 
eigen  value  problem.  The  configuration  of  the  airfoil  pitching  at  a  constant  rate  about 
its  quarter  chord  point  is  shown  in  figure  82.  X-Y  is  the  reference  frame  attached  to  the 
airfoil  at  the  quarter  chord  point  and  aligned  along  the  chord.  Xi-Yi  is  the  reference 
frame  attached  to  the  airfoil  surface  at  {xo,  Vo)  and  aligned  along  the  surface.  Boundary 
conditions  need  to  be  determined  at  points  A  and  B  on  the  Yi  axis.  A  is  a  point  on  the 
airfoil  surface,  while  B  is  a  point  on  the  Yi  axis  at  a  large  distance  from  the  airfoil.  The 
no-slip  condition  can  be  applied  at  the  ‘A’  boundary,  which  is  written  as 

ui  =  vi  =  0  at  li  =  0  (200) 

ui  and  vx  are  the  velocities  in  the  Xx-Yx  frame  of  reference.  Therefore  the  boundary 
condition  applied  at  ‘A’  is 

^p{yx)  =  D(p{yx)  =  0  at  Yi  =  0  (201) 

To  determine  the  boundary  condition  at  point  ‘B’'*  on  the  Yx  axis,  let  us  consider  the 
velocities  in  the  Xx-Yx  frame  at  that  point 

ux  =  Uoo  cos  6  —  ily  cos  9  -f  Q,x  sin  6  (202) 

ui  =  —Uoo  sin  0  fly  sin  ^  4-  fla:  cos  6  (203) 

corresponds  to  a  point  where  14  =  oo 
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Figure  82:  Configuration  of  the  airfoil 


But 

S/  =  yo  +  2/1  cos  6  and  x  =  Xo  —  yi  sin  9 
Therefore  equations  202  and  204  become 

«i  =  Uoo  cos  9  —  Q,  cos  9{yo  +  t/i  cos  ^)  +  sin  9{xo  —  y\  sin  9) 
=  Uoo  cos  9  —  Uyo  cos  9  +  Qxo  sin  9  —  fiyi 


(204) 


(205) 


vi  = —Uoo  sin  0  sin  0(^0  +  2/i  cos  ^)  +  cos  ^(aio  —  2/1  sin  ^) 
=  —Uoo  sin  9  +  Uyo  sin  9  +  Uxo  cos  9 

Therefore 

d^Ui  dv\ 

~JT  =  ^ 

dyi  dyi 

Dip{yi)  =  D^ip{yi)  =  0  at  Yi  =  oo 

To  summarize,  the  boundary  conditions  applied  at  the  two  boundaries  are 

(p{yi)  =  D(p{yi)  =0  at  yi  =  0  and 

D(p{yi)  =  D^(p{yi)  =  0  at  yi  =  oo 


(206) 

(207) 

(208) 

(209) 

(210) 
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F  Results  using  Implicit  Unstructured 
Navier- Stokes  Code 

During  the  period  of  the  no-cost  extension  (30  March  1995  -  30  March  1996),  the 
research  effort  focused  on  the  development  and  implementation  of  a  preconditioner  for 
the  implicit  unstructured  Navier-Stokes  algorithm,  and  the  application  of  the  algorithm 
to  the  simulation  of  flow  past  a  pitching  airfoil.  Two  cases  were  considered.  The  first 
case  {Rtc  —  10^,  M^o  —  0.2,  =  0.2)  is  identical  to  Case  1  (see  Sections  9.1  and  9.2), 

and  thus  provides  a  validation  of  the  implicit  unstructured  Navier-Stokes  algorithm. 
The  second  case  [Re^  =  2  x  10^,  M^o  =  0.2,  Q,^  =  0.2)  represents  a  configuration  not 
previously  simulated,  and  provides  new  insight  into  the  effects  of  Reynolds  number  on 
incipient  separation  for  pitching  airfoils.  Specifically,  the  primary  recirculation  region 
appears  at  a  =  14.5°  for  Rcc  =  2  x  10^,  compared  to  a  =  14.99°  for  Rbc  =  10^. 
Furthermore,  the  primary  recirculation  region  appears  closer  to  the  leading  edge  for  Rcc  = 
2  X  10^.  These  effects  of  increasing  Rcc  are  similar  to  those  obtained  at  Moo  =  0.5  for 
Rcc  =  10^  and  10®  at  =  0.2  (z.e..  Cases  4  and  7  in  Sections  9.5  and  9.8,  respectively). 


F.l  Numerical  Algorithm 

The  implicit  algorithm  for  the  Navier-Stokes  equations  using  an  unstructured  grid 
of  triangles  is  described  in  Section  5.  The  method  is  second-order  accurate  in  space, 
and  either  first-order  or  second-order  accurate  in  time,  depending  on  the  choice  of  the 
temporal  integration  parameters.  For  these  computations,  the  first-order  time  accurate 
Euler  method  is  used  with  the  time  step  At  chosen  to  be  the  same  as  employed  for  the 
structured  grid  computation  (see  Section  9.1).  Additionally,  the  Jacobian  is  updated 
every  20  time  steps  to  reduce  the  CPU  time.  Thus,  the  overall  algorithm  is  first-order 
accurate  in  time.  Nonetheless,  the  results  described  herein  indicate  that  a  highly  accurate 
solution  is  achieved,  despite  the  use  of  the  first-order  time  accuracy,  due  to  the  small  value 
of  the  time  step  At. 
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Preconditioning  was  developed  for,  and  implemented  in,  the  BiCGSTAB  algorithm 
(Section  5.7)  using  an  incomplete  LU  factorization  of  the  Jacobian  matrix  as  the  pre¬ 
conditioner.  An  LU  factorization  of  the  Jacobian  is  performed,  but  at  each  stage  of  the 
factorization  the  fiU  elements  are  ignored.  Thus  the  incomplete  LU  factorization,  to  be 
used  as  the  preconditioner,  has  the  same  sparsity  structure  as  the  original  Jacobian  ma¬ 
trix.  This  incomplete  factorization  is  known  as  the  ILU(O)  factorization,  since  no  fiU-in 
is  allowed.  The  preconditioning  matrix  in  its  LU  form  is  then  used  to  solve  systems  of 
equations  as  required  in  BiCGSTAB. 

F.2  Details  of  Computations 

An  unstructured  grid  based  on  13319  points  used  by  Choudhuri  [4],  and  having  26278 
triangular  cells  was  used.  There  were  180  points  on  the  airfoil  surfaces.  The  grid  for 
Rtc  =  2  X  10^  was  derived  from  the  grid  for  Rcc  —  10'*  by  moving  the  points  1  / \/2  closer 
to  the  airfoil.  Details  of  the  computational  grid  can  be  found  in  Table  1  and  Figs.  1  and 
2.  All  distances  are  normalized  by  the  airfoil  chord  c.  The  grid  is  attached  to  the  airfoil. 


Table  1:  Details  of  Computational  Grid 


Reynolds  Number 

Number  of  cells 

26278 

26278 

Minimum  As 

Maximum  As 

Minimum  An 

Maximum  An 

Nbl 

10 

10 

As  Distance  along  airfoil 

An  Distance  normal  to  airfoil 

Nbl  Number  of  points  in  boundary 


layer  at  mid-chord 

Amax,  ^min  Maximum,  minimum  area 

Results  were  obtained  for  laminar  flow  over  a  pitching  NACA  0012  airfoil  at  a  Mach 
number  of  0.2,  pitching  rate  fip  of  0.2  and  Reynolds  numbers  of  10*  and  2  x  10*.  The 
Jacobian  was  updated  every  20  timesteps,  with  a  timestep  of  5  x  10“®.  Computations 
were  performed  on  a  CRAY  C90  supercomputer.  As  described  previously,  the  grid  is 
attached  to  the  airfoil,  and  therefore  rotates  with  the  airfoil.  The  velocity  is  computed 
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in  the  inertial  frame;  for  analysis,  the  rotational  motion  of  the  airfoil  is  removed  and 
the  velocities  computed  in  the  airfoil  frame.  The  streamhnes  are  then  plotted  using  the 
velocity  components  in  the  airfoil  frame. 

F.3  Results  for  Rcc  =  10'^ 

The  Rec  =  10^  case  is  identical  to  Case  1  (see  Sections  9.1  and  9.2),  and  provides  a 
validation  of  the  imphcit  unstructured  grid  algorithm  for  unsteady  flows. 

The  computed  Hft,  drag  and  moment  coefiicients,  shown  in  Fig.  3  are  in  close  agree¬ 
ment  with  the  previous  computations  using  both  the  exphcit  unstructured  grid  and  im¬ 
phcit  structured  grid  Navier-Stokes  algorithms  (see  [4],  Fig.  3). 

In  Figs.  4  to  11,  13  and  14,  the  instantaneous  streamhne  plots  are  shown  for  both  the 
exphcit  unstructured  grid  and  imphcit  unstructured  grid  computations  at  14.5°,  16.5°, 
19.5°,  21.5°  and  22.5°.  The  imphcit  and  exphcit  algorithms  are  in  excellent  agreement. 
It  was  previously  shown  [4]  that  the  exphcit  unstructured  and  imphcit  structured  grid 
algorithms  displayed  excellent  agreement  for  this  case. 

The  algorithm  required  14.6  hours  of  CPU  time  on  the  CRAY  C90  to  pitch  the  airfoil 
up  to  22.5  degrees.  The  CPU  time  per  timestep  was  about  350  seconds  for  a  step  with 
an  update  of  the  Jacobian  and  the  preconditioner,  and  7  seconds  for  a  step  without  the 
update.  Five  iterations  were  required  for  convergence  of  BiCGSTAB  at  each  timestep. 

F.4  Results  for  Rcc  =  2  x  10'* 

The  Rcc  =  2  X  10^  represents  a  new  case.  The  flow  development  is  similar  to 
Rcc  =  10^,  with  a  few  major  differences.  As  shown  in  Fig.  12,  the  hft  and  drag 

coefficient  generally  increase  with  the  angle  of  attack;  the  moment  coefficient  is  some¬ 
what  constant.  Compared  to  Rcc  =  10^,  the  hft  is  shghtly  increased  and  drag  is  shghtly 
decreased.  The  variation  of  hft  and  drag  coefficients  with  angle  of  attack  is  not  as  mono¬ 
tonic  at  the  higher  Reynolds  number.  These  differences  are  illustrated  in  Fig.  15. 

The  flow  development  is  quahtatively  similar  to  Rtc  =  10^:  the  appearance  of  the  pri¬ 
mary  recirculating  region  at  14.5°  (Fig.  16),  its  growth  at  16.5°  (Fig.  17),  the  appearance 
of  a  secondary  recirculating  region  by  19.5°  (Fig.  18),  the  formation  of  a  tertiary  recir¬ 
culating  region  by  22.5°  (Fig.  19),  and  the  appearance  of  a  fourth  recirculating  region 
by  23.5°  (Fig.  20). 
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An  examination  of  the  angles  at  which  the  primary  recirculating  region  forms  shows 
that  increasing  Rtc  accelerates  the  onset  of  flow  separation,  with  the  recirculation  region 
appearing  at  a  lower  angle  of  attack  for  the  higher  Rtc-  The  recirculation  region  also 
appears  to  be  closer  to  the  leading  edge.  These  trends  are  in  agreement  with  pubhshed 
results  [3]  at  M^o  =  0.5  and  =  0.2.  At  this  higher  Mach  number,  an  increase  in  Rcc 
from  10^  to  10®  accelerates  the  formation  of  the  primary  recirculation  region  from  18.8° 
to  14.9°  and  moves  it  forward. 

The  algorithm  required  15.7  hours  of  CPU  time  on  the  CRAY  COO  to  pitch  the  airfoil 
up  to  22.5  degrees.  The  CPU  time  per  timestep  was  about  350  seconds  for  a  step  with 
an  update  of  the  Jacobian  and  the  preconditioner,  and  8  seconds  for  a  step  without  the 
update.  At  each  timestep,  six  BiCGSTAB  iterations  were  required  for  convergence  of  the 
iterative  solution  of  the  matrix  equations. 
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Figure  1:  NACA0012  Airfoil  Mesh  Rcc  —  10“* 


Figure  2:  NACA0012  Airfoil  Mesh  Rcc  =  2  x  10^ 


Figure  3:  NACA0012  Pitching  Airfoil  Rcc  =  10“^,  Lift, Drag  and  Moment  Coefficients 
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Figure  4:  NACA0012  Pitching  Air-  Figure  7:  NACA0012  Pitching  Air¬ 
foil,  Implicit  Computation,  Rcc  =  10"* ,  foil,  Explicit  Computation,  Rbc  =  10^, 

Streamlines,  a  =  14.5°  Streamlines,  a  =  14.5° 


Figure  5:  NACA0012  Pitching  Air-  Figure  8:  NACA0012  Pitching  Air¬ 
foil,  Implicit  Computation,  Rcc  =  10^,  foil,  Explicit  Computation,  Rsc  =  10^, 

Streamlines,  a  =  16.5°  Streamlines,  a  =  16.5° 


Figure  6:  NACA0012  Pitching  Air-  Figure  9:  NACA0012  Pitching  Air¬ 
foil,  Implicit  Computation,  Rcc  =  10^,  foil,  Explicit  Computation,  Rcc  =  10"*, 

Streamlines,  a  =  19.5°  Streamlines,  a  =  19.5° 
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Figure  10:  NACA0012  Pitching  Air¬ 
foil,  Implicit  Computation,  Rtc  =  10'*, 
Streamlines,  a  =  21.5° 


Figure  13:  NACA0012  Pitching  Air¬ 
foil,  Explicit  Computation,  Rcc  =  10“*, 
Streamlines,  a  =  21.5° 


Figure  11:  NACA0012  Pitching  Air¬ 
foil,  Implicit  Computation,  Rtc  =  10^, 
Streamlines,  a  =  22.5° 


Figure  14:  NACA0012  Pitching  Air¬ 
foil,  Explicit  Computation,  Rcc  =  10'^, 
Streamlines,  a  =  22.5° 


Figure  12:  NACA0012  Pitching  Airfoil, 
Explicit  Computation,  Rcc  —  2  x  10^, 
Lift,  Drag  and  Moment  Coefficients 
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Figure  15:  NACA0012  Pitching  Airfoil 
Comparison  of  Lift,  Drag  and  Moment 
Coefficients  at  Rcc  =  10^  and  2  x  10^ 


=  14.5“ 


Figure  16:  NACA0012  Pitching  Airfoil, 
Implicit  Computation,  Rcc  =  2  x  10^, 
Streamlines,  a  =  14.5° 


Figure  17:  NACA0012  Pitching  Airfoil, 
Implicit  Computation,  Rec  =  2  x  10^, 
Streamlines,  a  =  16.5° 


Figure  19:  NACA0012  Pitching  Airfoil, 
Implicit  Computation,  Re^  =  2  x  10^, 
Streamlines,  a  —  22.5° 


Figure  20:  NACA0012  Pitching  Airfoil, 
Implicit  Computation,  Re^.  =  2  x  10^, 
Streamlines,  a  =  23.5° 


Figure  18:  NACA0012  Pitching  Airfoil, 
Implicit  Computation,  Rcc  =  2  x  10^, 
Streamlines,  a  =  19.5° 


142 


